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ABSTRACT 

'  The  offshore  tide  becomes  distorted  as  it  propagates  into  shallow  inlet/ 
estuarine  systems.  Time  asymmetries  develop  in  the  rise  and  fall  of  sea 
surface  with  consequent  time  and  magnitude  asymmetries  in  tidal  currents. 
Flood-dominant  estuaries  are  characterized  by  longer  falling  tides  and 
stronger  flood  currents  while  ebb-dominant  estuaries  have  longer  rising 
tides  and  stronger  ebb  currents.  The  asymmetries  are  reflected  in  the  non¬ 
linear  growth  of  harmonics  and  compound  tides  of  the  principal  equilibrium 
tidal  constituents.  This  dissertation  consists  of  three  papers  which 
examine  the  development  of  tidal  asymmetries  in  shallow  estuarine  systems: 
a  study  of  the  recent  migration  history  of  Nauset  Inlet  (MA),  a  shallow 
estuarine  system  located  on  Cape  Cod;  an  analysis  of  the  results  of  a  series 
of  field  experiments  conducted  at  Nauset;  a  numerical  model  study  of  the 
types  of  estuarine  characteristics  controlling  tidal  asymmetry.  The  ana¬ 
lysis  of  field  results  focuses  on  sea  surface  measurements.  Non-linear 
distortion  of  the  tide  at  Nauset  is  characterized  by  the  strong  growth  of 
harmonics  and  compound  constituents  particularly  in  the  quarter-diurnal 
band.  Phase  relationships  between  the  forced  constituents  and  their  parents 
produce  a  flood-dominant  estuary.  Numerical  modeling  of  M2  tidal  propaga¬ 
tion  in  shallow  estuarine  channels  utilizes  the  one-dimensional  equations 
of  motion. y Shallow,  frictionally  dominated  channels  with  moderate  tidal 
flat  area  develop  a  flood-dominant  asymmetry  while  deeper  channels  with 
extensive  tidal  flats  develop  an  ebb-dominant  asymmetry.  Model  results 
are  supported  by  observations  of  tidal  asymmetry  in  natural  estuaries. 
Implications  of  non-linear  tidal  distortion  on  bedload  and  suspended 
material  transport  are  profound.  Flood-dominant  estuaries  tend  to  import 
sediment  if  the  supply  is  adequate  whereas  ebb-dominant  estuaries  can  flush 
entering  sediment  effectively.  Over  long  time  periods,  flood-dominant 
estuaries  may  eventually  fill.  Ebb-dominant  estuaries  may  represent  more 
stable  long-term  configurations. 

'  *  •  - 
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CHAPTER  I 
INTRODUCTION 
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Shallow  estuarine/tidal  marsh  systems  connected  to  the  ocean  by  a 
narrow  inlet  are  ubiquitous  features  along  barrier  beach  coastlines. 

Barrier  beaches  and  their  associated  tidal  inlets  comprise  approximately 
13%  of  the  world's  coastline  (King,  1972),  most  occurring  in  areas  with 
low  to  moderate  tidal  range.  The  longest  stretch  of  barrier  beaches  and 
tidal  inlets  in  the  world  is  located  along  the  U.S.  East  and  Gulf  of  Mexico 
coasts,  affording  protection  to  over  5000  km  of  navigable  waters.  Inlet/ 
estuarine  systems  along  these  coastlines  play  important  roles  in  the  bio¬ 
logical  cycles  of  many  organisms  and  in  nutrient  exchange  with  coastal 
waters.  In  addition,  many  are  used,  or  are  intended  for  use,  as  naviga¬ 
tional  channels.  As  the  offshore  tide  propagates  into  these  systems,  it 
can  become  distorted  as  a  result  of  non-linear  fluid  mechanical  proces¬ 
ses.  The  tidal  distortion  takes  the  form  of  time  asymmetries  in  the  rise 
and  fall  of  sea  surface  and  consequent  time /magnitude  asymmetries  in  tidal 
currents.  Patterns  of  sediment  transport  and  nutrient/pollutant  exchange 
can  be  strongly  affected  by  the  tidal  distortion.  The  long  term  stability 
and  evolution  of  these  estuaries  are  partly  controlled  by  their  response 
to  tidal  forcing. 

Coastal  geologists  have  focused  on  the  velocity  asymmetries  frequently 
observed  in  the  tidal  inlet  channel  and  classified  inlets  as  either  "ebb 
dominant"  or  "flood  dominant"  (eg.,  Hayes,  1975).  Many  inlet  channels  are 
termed  "ebb  dominant"  because  it  is  observed  that  ebb  currents  continue  to 
flow  even  after  the  ocean  tide  has  turned  to  rise.  As  a  result,  flood  cur¬ 
rents  are  confined  to  shallow  channels  on  either  side  of  the  main  inlet 
channel  during  early  stages  of  rising  tide  (see  figure  1  for  inlet  morpho¬ 
logy).  This  situation  is  caused  by  lags  in  low  water  from  the  ocean  to  the 
far  reaches  of  the  estuary.  The  terms  "ebb  dominant"  and  "flood  dominant" 


EBB  TIDAL  DELTA 


Figure  1.  Sketch  of  dominant  morphologic  features  in  a  tidal  inlet. 
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have  generally  been  applied  to  currents  at  the  inlet  as  opposed  to  further 
within  the  estuary.  Although  most  work  has  focused  on  the  time-velocity 
asymmetry  of  tidal  currents  investigators  have  noted  the  time  asymmetry 
present  in  sea  surface  (eg.,  Byrne  et  al.,  1975;  Boon  and  Byrne,  1981; 
FitzGerald  and  Nummedal,  1983).  The  duration  asymmetry  in  rising  and 
falling  tide  can  impart  an  ebb  or  flood  dominance  to  the  estuary  as  a 
whole.  For  example,  an  estuary  characterized  by  a  longer  rising  tide  will 
tend  to  have  stronger  ebb  currents  ("ebb  dominant").  Estuaries  of  both 
types  occur  along  the  U.S.  East  coast  (Table  1). 

The  importance  of  tidal  asymmetries  to  processes  such  as  sediment 
transport  (eg.,  Postma,  1967;  Pingree  and  Griffiths,  1979)  has  long  been 
recognized.  The  presence  of  a  distorted  tidal  velocity  signal  can  produce 
a  net  suspended  and  bed  load  transport.  Within  estuaries,  net  patterns  of 
sediment  transport  can  have  important  implications  for  inlet/estuarine 
stability.  An  estuary  characterized  by  a  longer  falling  tide  and  hence 
stronger  flood  currents  may  show  a  pattern  of  net  up-estuary  sediment 
transport.  Such  an  estuary  may  be  unable  to  flush  effectively  sediments 
entering  the  system  through  the  inlet.  As  a  result,  the  estuary  could  be 
expected  to  gradually  fill,  with  consequent  reduction  of  tidal  prism  and 
development  of  channel  shoaling  problems.  The  tidal  response  of  these 
systems,  therefore,  may  play  an  important  role  in  their  geological  develop¬ 
ment.  Some  investigators,  in  fact,  have  hypothesized  that  an  inlet/estuar¬ 
ine  system  could  evolve  from  a  flood-dominant  transport  pattern  in  its 
early  stages  to  ebb-dominant  as  the  estuary  filled  (Byrne  and  Boon,  1976). 
Such  a  development  implies  a  feedback  between  estuarine  geometry  and  tidal 


fore ing . 
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Table  1 

Examples  of  Tidal  Asymmetry  along  U.S.  east  coast 


Longer  falling  tide  ("flood  dominant") 
Murrells  Inlet,  S.C. 

Nauset  Inlet,  MA 
Ogunquit  River  Inlet,  ME 


Reference 

Boon  and  Byrne,  1981 
Aubrey  and  Speer,  1983 
Fitzgerald  et  al.  (in  press) 


Longer  rising  tide  ("ebb  dominant") 
Wachapreague  Inlet,  VA 
North  Inlet,  S.C. 

Price  Inlet,  S.C. 


Boon  and  Byrne,  1981 
Nummedal  and  Humphries,  1978 
Fitzgerald  and  Nummedal,  1983 
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Over  shorter  time  scales,  many  investigators  have  been  interested  in 
patterns  of  net  suspended  material  transport.  This  is  particularly  impor¬ 
tant  for  biologists  concerned  whether  inlet/marsh  systems  expert  dissolved 
and  particulate  nutrients  into  coastal  waters  (Valiela  et^  al . ,  1978;  Teal, 
1983).  Nutrient  circulation  depends  in  part  on  transformations  which  occur 
in  the  tide  as  it  propa  as  into  these  shallow  systems.  Net  patterns  of 
suspended  material  transport  have  been  explained  by  means  of  the  time- 
velocity  asymmetries  developed  in  estuarine  channels  (eg.  Boon,  1975; 

Ward,  1981). 

Despite  the  large  number  of  problems  studied  in  shallow  estuaries,  few 
have  focused  on  the  tidal  propagation  problem  itself.  However,  it  is  the 
interaction  of  estuarine  geometry  and  tidal  forcing  which  produces  the 
observed  asymmetries  in  tidal  currents  and  the  rise  and  fall  of  sea 
surface.  A  number  of  investigators  have  utilized  numerical  models  to 
investigate  the  role  of  estuarine  hypsometry  in  distorting  the  tide  (Mota- 
Oliveira,  1970;  Seelig  and  Sorensen,  1978;  Boon  and  Byrne,  1981).  These 
studies  have  not  been  sufficiently  systematic.  In  addition,  they  did  not 
focus  on  the  spectral  response  of  the  model  estuaries  to  tidal  forcing. 

This  is  important  because  the  non-linear  process  of  tidal  distortion 
observed  in  the  field  reveals  itself  in  the  growth  of  harmonics  of  the 
principal  astronomic  constituents. 

Data  sets  obtained  in  the  field  have  largely  consisted  of  12  or  2b  hour 
cycles  of  velocity  measurements  typically  sampled  at  hourly  intervals. 

These  measurements  frequently  have  been  concentrated  near  the  inlet  channel 
(eg.,  Boothroyd  and  Hubbard,  1975;  Finley,  1975;  Hine,  1975;  Nummedal  and 
Humphries,  1978;  FitzGerald  and  Nummedal,  1983).  Long-term  records  of  sea 
surface  measurements  (sufficient  for  harmonic  analysis)  are  sparse,  with 
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only  a  few  comprehensive  surveys  (e.g.,  NOS  experiment,  Murrells  Inlet, 

SC,  Boon  and  Byrne,  1981).  A  comprehensive  field  program  supported  by 
modeling  is  required  to  understand  tidal  propagation  in  these  shallow 
systems.  The  key  questions  are:  What  are  the  kinematics  and  mechanics  of 
the  non-linear  distortion  of  the  tide  within  a  shallow  estuary?  What  are 
the  geometry  and  friction  characteristics  of  estuaries  which  produce  the 
different  types  of  asymmetry?  The  present  study  addresses  these  questions 
It  has  involved  a  field  study  at  Nauset  Inlet,  MA,  a  tidal  inlet/marsh 
complex  representative  of  many  such  features  along  the  U.S.  east  coast. 

In  addition,  numerical  modeling  was  employed  as  a  diagnostic  tool  to 
investigate  the  important  aspects  of  tidal  propagation  identified  by  the 
field  experiments. 

The  thesis  consists  of  three  papers  (co-authored  with  D.  Aubrey):  a 
study  of  the  recent  migration  history  of  Nauset  Inlet;  an  analysis  of  the 
results  of  a  series  of  field  experiments  conducted  at  Nauset;  a  numerical 
model  study  of  the  types  of  estuarine  characteristics  responsible  for 
producing  different  tidal  asymmetries.  The  first  paper  places  the  field 
experiments  within  the  framework  of  the  recent  history  of  the  inlet/ 
estuarine  system.  The  Nauset  estuary/marsh  complex  is  serviced  by  an 
unstable  inlet  which  has  migrated  as  much  as  2  km  over  the  past  30  years. 
The  inlet  is  presently  moving  north  at  approximately  100  m/yr  and  is  more 
than  2  km  from  an  historically  preferred  location.  Such  large-scale  migra 
tion  changes  the  lengths  of  estuarine  tidal  channels  and  can  affect  its 
tidal  response.  This  paper  examines  the  migration  episodes  and  formulates 
three  hypotheses  explaining  the  observed  behavior. 

The  second  paper  discusses  results  from  field  experiments  designed  to 
examine  tidal  propagation  in  this  shallow  estuary.  The  analysis  focuses 
on  sea  surface  elevation  measurements,  supplementing  them  with  velocity 
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records  which  illustrate  the  type  of  tidal  asymmetry  developed  at  Nauset. 

The  non-linear  distortion  of  the  tide  is  characterized  by  the  growth  of 
harmonics  and  compound  constituents  of  the  primary  astronomic  frequencies. 

Sea  surface  records  of  sufficient  length  to  examine  this  problem  spectrally 
are  far  easier  to  obtain  than  velocity  records.  In  fact,  difficult  environ¬ 
mental  conditions  (shallow  water,  strong  flows,  extreme  sediment  transport 
rates)  preclude  long  velocity  records  from  most  regions  of  the  estuary. 

Hence  velocity  is  used  as  supplementary  information  to  a  detailed  investi¬ 
gation  of  sea  surface  records. 

The  numerical  section  presents  results  from  one-dimensional  (cross- 
sectionally  averaged)  modeling  of  estuarine  characteristics  contributing 
to  different  forms  of  tidal  asymmetry.  The  modeling  is  diagnostic  and 
does  not  attempt  to  predict  or  hindcast  conditions  at  Nauset.  These 
systems  are  physically  complex  and  require  two-dimensional  models  to 
hindcast  sea  surface  and  depth-mean  velocity.  Such  models  exist  (eg., 

Masch  et:  al.,  1977;  Butler,  1980)  but  they  do  not  explain  what  drives 
estuaries  to  different  tidal  responses.  In  addition,  their  ability  to 
represent  non-linear  estuarine  dynamics  accurately  is  unclear  since 
investigators  have  not  provided  spectral  comparisons  of  field  data  with 
model  solutions.  The  one-dimensional  models  are  used  to  examine  the  role 
of  channel  cross-section  shape,  tidal  flats  and  variable  levels  of  friction 
in  distorting  the  estuarine  tide.  Model  results  indicate  systems  with 
strongly  time-variable  channel  cross-section  and  moderate  extent  of  tidal 
flats  tend  to  be  flood-dominant.  Conversely,  estuaries  with  extensive 
tidal  flats  and  less  time-variability  in  channel  geometry  are  ebb-dominant. 
Results  of  the  field  experiments  at  Nauset  Inlet  and  observations  of  tidal 
asymmetry  in  other  shallow  estuaries  support  the  findings  of  the  model  study. 
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ABSTRACT 

Three  mechanisms  are  responsible  for  tidal  inlet  migration  in  an 
updrift  direction  (counter  to  net  longshore  transport):  1)  attachment  of 
distal  ebb  tide  delta  bars  to  the  downdrift  barrier  spit;  2)  storm-induced 
breaching  and  subsequent  stabilization  to  form  a  new  inlet;  and  3)  ebb  tide 
discharge  around  a  channel  bend  creating  a  three-dimensional  flow  pattern 
which  erodes  the  outer  channel  bank  and  accretes  on  the  inner  channel  bank. 
The  last  two  mechanisms  can  result  in  either  updrift  or  downdrift  inlet 
migration,  depending  on  channel  geometry  in  the  bay  and  barrier  beach  con¬ 
figuration.  The  last  mechanism  is  discussed  here  for  the  first  time. 
Analysis  of  historical  charts  and  aerial  photographs,  combined  with  an 
historical  storm  synthesis,  shows  that  all  three  mechanisms  are  active  at 
a  natural  tidal  inlet  along  a  sandy  coast  (Nauset  Inlet,  Cape  Cod,  MA). 

On  a  time  scale  of  ten  years,  these  mechanisms  were  effective  in  producing 
an  updrift  migration  of  more  than  two  kilometers.  Initiation  of  updrift 
migration  coincided  with  a  marked  increase  in  storm  frequency  perturbing 
the  historically  stable  inlet  position.  Subsequent  updrift  migration 
resulted  from  ebb-delta  bypassing  and  channel  bend  flows. 
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INTRODUCTION 

Migration  of  tidal  inlets  and  the  associated  changes  in  adjacent 
barrier  beaches  have  profound  implications  on  both  the  geological  evolution 
of  inlet/estuary  systems  and  the  short-term  stability  of  these  features. 

Past  studies  have  documented  many  instances  of  inlets  migrating  in  the 
direction  of  net  littoral  drift  along  sandy  shores,  but  have  uncovered  few 
cases  (e.g.,  Indian  River  Inlet,  Delaware,  and  Thorsminde  Inlet,  Denmark) 
where  inlets  appear  to  migrate  in  directions  opposed  to  the  dominant  long¬ 
shore  transport  direction  (Bruun,  1978).  Migration  of  tidal  inlets  in  any 
direction  accelerates  inlet-induced  changes  in  the  estuary.  The  estuary 
may  fill  in  with  littoral  sands  derived  from  updrift  sources,  as  flood  tide 
delta  growth  accompanies  the  migration  of  the  inlet,  and  marsh  development 
(colonization  and  plant  emergence)  becomes  more  variable  and  less  permanent. 

Previous  attempts  to  explain  a  reversal  in  direction  of  inlet  migration 
suggest  a  change  in  direction  of  net  littoral  drift,  causing  a  change  in 
migration  direction.  This  explanation  is  not  realistic  for  some  inlets 
where  wave  forcing  and  nearshore  bathymetry  have  remained  constant  through 
time.  This  study  presents  three  alternatives  to  explain  the  tendency  of 
some  inlets  to  migrate  updrift,  each  supported  by  historical  observations 
at  a  site  with  a  large-volume,  directionally-biased  littoral  drift. 

The  study  site  is  located  on  the  Atlantic  coast  of  Cape  Cod,  Massa¬ 
chusetts  (figure  1),  exposed  to  open  ocean  waves  from  the  east  and  a  two- 
meter  ocean  tide.  Offshore  bathymetry  and  sediments  are  described  else¬ 
where  (Aubrey,  Twichell  and  Pfirman,  1982).  Longshore  transport  rates  and 
directions  were  studied  by  Zeigler  (1954,  I960),  and  net  littoral  drift 
has  been  estimated  at  250,000  m3  per  year  towards  the  south  (U.S.  Army 
Corps  of  Engineers,  1969).  Sediment  is  derived  from  erosion  of  sea  cliffs 
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border  ing  Nauset  Inlet  to  the  north.  Overwash  processes  along  Nauset 
barrier  beaches  are  described  in  Zaremba  and  Leatherman  (1979).  Aubrey  and 
Speer  (1983,  1984)  discuss  tidal  flows  and  sediment  transport  in  the  bay 
and  inlet.  Other  sedimentologic  studies  of  the  region  are  found  in  a 
summary  volume  by  Leatherman  (1979). 

Sea  level  here  has  risen  an  average  of  about  3.5  mm/year.  This  rate 
is  three  times  greater  than  the  mean  sea  level  rise  of  one  mm/year  over 
the  past  2100  years  established  from  measurements  of  salt  marsh  peat 
accumulation  at  Barnstable  Harbor,  Cape  Cod  (Redfield  and  Rubin,  1962). 

As  discussed  by  Aubrey  and  Emery  (1983)  and  others,  although  short-term 
mean  sea  level  records  exhibit  considerable  oscillations  about  a  mean 
trend,  the  mean  rate  over  the  past  100  years  has  not  changed  significantly. 
Sea-level  rise  favors  landward  migration  of  evolving  barrier  beaches. 

METHODS 

Historical  charts  and  aerial  photographs  of  the  Nauset  Inlet  area 
dating  from  1670  and  1938,  respectively,  were  examined  to  define  and  (where 
possible)  quantify  changes  in  inlet  position  and  morphology.  Historical 
data  (figure  2)  were  obtained  from  a  variety  of  sources  including  govern¬ 
ment  agencies,  the  National  Archives,  the  Library  of  Congress,  the  Woods 
Hole  Oceanographic  Institution,  and  private  industry  (Speer,  Aubrey  and 
Ruder,  1982;  Appendices  1  and  2).  Chart  coverage  is  dense  from  1790  to 
present  (coverage  was  sparse  before  1790),  and  good  aerial  photographic 
coverage  exists  from  1951  to  present  (only  one  aerial  photo  sequence  was 
available  prior  to  1951,  taken  in  1938).  Small  scale  and  uncertain  mapping 
techniques  used  in  pre-1846  historical  charts  make  it  difficult  to  quantify 
changes  in  inlet  morphology  during  this  period,  but  these  charts  were  valu¬ 
able  in  depicting  general  trends  in  inlet  morphology.  Care  was  required  in 
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Figure  2.  Summary  of  Nauset  Inlet  aerial  photography  (125  sets  since  1938) 
and  historical  charts  (120  sets  since  1670). 
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interpreting  the  charts  because  several  of  the  charts  from  the  1800 's  did 
not  specify  survey  dates,  and  were  merely  reproductions  of  earlier  and 
perhaps  outdated  surveys.  Also,  in  the  case  of  U.S.  Coast  and  Geodetic 
Survey  (USC&GS)  charts,  only  limited  shoreline  segments  were  updated 
between  editions. 

Aerial  photographs  provide  more  detailed  information  than  the  charts 
because  they  are  generally  larger  in  scale  (allowing  resolution  of  shore¬ 
line  features  such  as  bars  and  marshes).  They  also  provide  more  comprehen¬ 
sive  temporal  coverage  for  a  limited  period  (1951  to  1981)  than  do  the 
charts,  and  the  dates  of  coverage  are  unambiguous.  Fifty  vertical  sets  of 
the  125  photographs  available  were  measured  to  quantify  inlet  and  spit 
migration  at  Nauset.  The  remaining  photographs  were  not  measured  because 
they  were  taken  at  oblique  angles,  were  poorly  fitted  mosaic  series,  or 
lacked  sufficient  ground  control  to  assure  measurement  accuracy.  However, 
they  were  instrumental  in  providing  a  continuous  record  of  relative  changes 
in  inlet  and  spit  locations  during  the  past  30  years. 

Measurements  of  spit  and  inlet  locations  are  relative  to  a  baseline, 
sub-parallel  to  the  shoreline,  established  between  well-defined,  permanent 
features  identified  on  each  set  of  aerial  photographs  (figure  3).  The 
known  length  of  this  baseline  provided  a  consistent  determination  of  scale 
for  all  photos.  Uncertainty  in  some  measurements  resulted  when  one  of  the 
two  primary  reference  points  was  absent  from  a  particular  photo  mosaic. 

In  these  cases,  secondary  landmarks  were  used  along  with  geometrical  rela¬ 
tions  to  define  the  baseline  from  the  one  available  endpoint.  As  a  result 
of  such  variations  in  the  photographs,  overall  accuracy  of  measurements  is 
estimated  to  be  +  15  m,  despite  a  measurement  resolution  of  5  m. 
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Figure  3.  Baseline  location  for  inlet  migration  measurements.  The  northern 

reference  point  adjacent  to  Coast  Guard  beach  is  Coast  Guard  Station. 
Distance  between  adjacent  reference  lines  (1-10)  is  535  m. 
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RESULTS 

Analysis  of  historical  charts  and  photographs  reveals  patterns  of  inlet 
migration  and  barrier  beach  elongation/shortening  at  Nauset.  Nauset  Inlet 
has  migrated  extensively  over  the  past  30  years,  related  in  part  to  an 
increase  in  storm  frequency,  as  discussed  below. 

INLET /BARRIER  BEACH  MIGRATION 

Historical  charts  (dating  from  1779)  and  aerial  photography  (dating 
from  1938)  show  the  preferred  inlet  location  to  have  been  just  north  of 
Nauset  Heights,  at  the  southern  (downdrift)  extremity  of  the  bay  drainage 
system  (figure  4).  None  of  the  charts  (up  to  1946)  depict  a  significant 
south  spit.  Although  aperiodic  historical  coverage  might  have  undersampled 
previous  episodes  of  inlet  migration,  the  persistence  of  a  southern  loca¬ 
tion  suggests  this  was  an  historically  stable  inlet  configuration. 

Aerial  photographs  from  1938  and  a  1946  USGS  chart  confirm  an  inlet 
location  just  north  of  Nauset  Heights,  with  no  south  barrier  apparent 
(figures  5  and  6).  From  the  1950's  into  the  early  1980's,  the  inlet  has 
been  active  with  three  distinct  cycles  of  northward  (updrift)  movement. 

The  first  two  of  these  (figure  5;  1952-1957;  1965-1972)  resulted  in  a  pat¬ 
tern  of  overlapping  spits.  In  both  cases,  the  length  of  the  north  spit 
remained  approximately  stable  while  the  south  spit  continually  grew  north. 
The  third  cycle  (figure  5;  1972-1984)  has  involved  substantial  erosion  of 
the  north  spit  along  with  northward  growth  of  the  south  sp't. 

The  first  cycle  of  northward  growth  was  initiated  by  storm  activity. 

The  north  barrier  was  breached  and  a  remnant  of  the  barrier  located  south 
of  the  new  breach  attached  to  the  south  barrier.  Subsequent  southward 
growth  of  the  north  barrier  through  attachment  to  an  island  in  the  marsh 
(Inlet  Marsh),  and  northward  growth  of  the  south  spit,  resulted  in  the  over- 
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baseline  from  Coast  Guard  Station  (see  figure  3  caption;,  btip] 
patterns  indicate  periods  when  South  Spit  overlapped  North  Spit 
and  extended  farther  north.  Vertical  lines  indicate  periods  fo 
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7.  Tracings  of  vertical  aerial  photographs  of  Nauset  Inlet  between 
1965  and  1978. 
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lapping  pattern  of  the  mid-1950's  (figure  6).  During  this  period  of 
barrier  growth,  the  elongated  inlet  channel  extended  to  the  north  from  its 
original  southern  location.  A  series  of  storms  in  the  late  1950's  and 
early  1960's  re-established  the  inlet  to  its  southernmost  position  immedi¬ 
ately  adjacent  to  Nauset  Heights  (figure  6f,  7a).  By  April  1965,  the  inlet/ 
barrier  beach  configuration  was  similar  to  that  in  1938,  with  the  sole  major 
difference  being  landward  migration  of  the  north  spit  due  to  repeated  storm 
overwashes  (Zaremba  and  Leatherman,  1979).  The  second  cycle  of  northward 
migration  was  characterized  by  extension  of  the  south  barrier  while  the 
north  barrier  remained  in  approximately  the  same  location.  Between  1965 
and  1972,  the  south  spit  extended  nearly  900  m  to  the  north  (figure  7b,  7c). 
As  in  the  mid-1950's,  the  inlet  channel  also  extended  north,  accompanying 
this  barrier  growth. 

The  final  cycle  of  northward  movement  was  initiated  by  a  storm  breach 
in  the  north  spit  in  the  spring  of  1972  (figure  7c).  This  resulted  in  the 
present  phase  of  northward  migration  which  is  qualitatively  different  from 
the  previous  two  episodes.  In  this  third  instance,  the  main  inlet  channel 
stabilized  in  the  location  of  the  newly  formed  breach  with  the  former  inlet 
closing  off.  Since  1972,  the  inlet  channel  and  the  south  barrier  have  been 
growing  northward  at  a  rate  of  approximately  100  m/yr.  Unlike  previous 
cycles,  the  north  spit  is  steadily  eroding  and  the  main  inlet  channel  is 
actually  moving  north,  instead  of  simply  extending  to  the  north  from  a 
southern  location  (figure  7d-7f).  No  major  breaches  affecting  the  stabil¬ 
ity  of  the  south  spit  have  occurred  during  this  latest  period.  A  large 
overwash  occurred  at  the  northernmost  part  of  the  north  spit  during  the 
6  February  1978  blizzard.  Since  this  overwash  emptied  into  a  shallow 
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(<1  m  deep),  broad  bay  (Nauset  Bay),  the  overwash  did  not  evolve  into  a 
permanent  breach.  Such  an  overwash  occurring  on  the  south  spit  would 
probably  result  in  a  new  inlet  position. 

STORM  ANALYSIS 

The  Atlantic  shore  of  Cape  Cod  is  frequently  buffeted  by  storms  which 
have  the  potential  to  cause  dramatic  changes  in  shoreline  configuration. 

A  U.S.  Army  Corps  of  Engineers  report  (1979)  cites  160  gales  with  wind 
speeds  greater  than  32  mph  between  1870-1975.  Half  of  these  were  north¬ 
easters.  Both  tropical  and  extra-tropical  (including  northeasters) 
cyclones  produce  dramatic  changes  at  Nauset  Inlet  because  of  the  geograph¬ 
ical  orientation  of  the  outer  Cape. 

Three  types  of  storm  data  were  collected  for  comparison  with  large- 
scale  morphologic  changes  at  Nauset  Inlet: 

a)  Hayden  and  Smith  (1982)  compiled  a  monthly  list  of  cyclone  occur¬ 
rences  off  the  east  coast  between  1885  and  1982,  using  as  a  data  base  the 
"Tracks  of  the  Centers  of  Cyclones  at  Sea  Level"  published  by  Monthly 
Weather  Review  and  in  recent  years  by  Mariners  Weather  Log.  Cyclone 
statistics  (both  tropical  and  extra-tropical)  are  available  on  2.5°  lati¬ 
tude  by  5°  longitude  grid  cells.  The  four  grid  cells  bordering  the  Cape 
Cod  region  to  the  east  and  southeast  (total  area  covered  is  60°W  to  70°W, 
37t°N  to  42t°N)  are  used  as  the  region  of  storm  influence  for  the 
study  area.  For  generation  of  year-by-year  and  monthly  mean  statistics, 
storm  values  for  the  four  grid  cells  are  summed.  Although  this  yields  an 
overestimate  of  the  number  of  storms  (the  same  storm  may  pass  more  than 
one  grid  cell),  it  will  still  provide  a  qualitative  indication  of  storm 
duration  and  persistence,  since  on  the  average  a  storm  tracking  through 
two  grid  cells  generates  waves  in  the  study  area  for  a  longer  period  of 
time  than  one  passing  a  single  grid  cell. 
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Cyclone  statistics  resulting  from  the  averaging  serve  as  a  crude  indi¬ 
cator  of  wave  activity.  Large  cyclone  counts  suggest  high  wave  activity; 
a  small  cyclone  count  represents  low  wave  activity.  Persistence  and 
frequency  of  storms  are  our  criteria  for  wave  intensity.  Clearly,  storm 
intensity  or  magnitude  would  be  a  useful  weighting  factor  for  linking 
waves  and  storms;  unfortunately,  this  information  is  not  available. 

The  period  from  1885  through  1949  experienced  a  relatively  low  incid¬ 
ence  of  storm  activity  (figure  8).  Within  this  low  background  level,  the 
periods  from  1885-1893,  1921-1924,  and  1930-1941  have  local  maxima  in 
cyclone  frequency.  The  last  thirty  years  of  the  record  show  consistently 
higher  cyclone  frequency,  with  local  maxima  at  1950-1954,  1961-1962,  1972 
and  1974.  Although  the  absolute  number  of  storms  may  be  sensitive  to  the 
quality  and  quantity  of  weather  observation  stations,  local  trends  (minima 
and  maxima)  are  valid  indicators  of  relative  storm  occurrence. 

b)  Another  source  of  storm  incidence  data  was  the  U.S.  Army  Waterways 
Experiment  Station  (WES)  wave  hindcast  program  (data  provided  by  W.  Birke- 
meier).  This  program  computes  nearshore  wave  height  statistics  based  on 
weather  observations  and  local  bathymetry.  The  study  identified  the  157 
largest  storm  events  from  1956  to  1976  (inclusive).  These  storms  were 
assigned  recurrence  intervals  according  to  their  rankings,  allowing  for 
weighting  of  storms  by  severity.  The  WES  compilation  (figure  9)  correlates 
well  with  cyclone  data.  W.E.S.  data  show  high  storm  activity  in  1956, 

1962,  and  1972;  however,  it  also  indicates  a  high  level  of  storm  activity 
in  1969  which  does  not  appear  in  cyclone  data.  Differences  between  the 
two  data  sets  are  the  result  both  of  weighting  procedures  and  different 
representations  of  the  data  base. 
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YEARLY  STORM  RECURRENCE  INTERVALS 

1956  -  1975 


Figure  9.  Yearly  storm  activity  from  1956-1976  off  Cape  Cod,  Massachusetts, 
determined  by  W.E.S.  hindcast. 
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Figure  10.  Compilation  of  storm  events  on  an  annual  basis  from  1800  through  1980. 

Sources  are  newspapers,  historical  descriptions,  and  published  tropical 
storm  tracks. 
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c)  Finally,  a  list  of  major  storms  affecting  the  outer  Cape  was  com¬ 
piled  from  newspapers,  historical  descriptions,  and  published  tropical 
storm  tracks  (figure  10).  This  list  is  incomplete  since  prior  to  1948  it 
only  includes  hurricanes  and  storms  of  historical  significance.  It  is 
possible  to  identify  specific  storms  which  are  likely  to  cause  changes  at 
Nauset  Inlet,  although  irregular  sampling  afforded  by  aerial  photography 
(figure  2)  makes  direct  correlation  difficult.  Through  this  method,  ten 
significant  storms  were  found  that  were  not  hindcast  in  the  WES  study. 

DISCUSSION 

Aerial  photographs  and  historical  charts  reveal  patterns  of  inlet/ 
barrier  beach  change  with  which  conceptual  models  must  be  consistent.  The 
important  features  of  Nauset  Inlet's  migration  patterns  are:  the  histor¬ 
ical  stability  of  the  southernmost  inlet  entrance;  the  role  of  storms  in 
initiating  major  changes  in  the  inlet/barrier  beach  system;  and  the  recent 
tendency  for  the  inlet  to  move  in  a  direction  opposite  the  predominant 
longshore  drift.  Migration  of  the  inlet  with  accompanying  changes  in  the 
barrier  beaches  takes  place  on  essentially  two  different  time  scales. 

Major  relocations  of  the  inlet,  involving  longshore  movements  of  hundreds 
of  meters  in  several  days,  occur  episodically  during  large  storms,  and  have 
a  recurrence  on  the  order  of  a  decade.  The  other  important  time  scale  is 
associated  with  the  recent  steady  migration  of  the  inlet  in  a  general  north¬ 
ward  direction,  responding  to  the  combined  effects  of  wave  activity,  tidal 
flows,  and  longshore  sand  transport.  The  magnitude  of  this  movement  is  on 
the  order  of  100  m/yr.  Northward  migration  of  the  inlet  is  accompanied  by 
extension  of  the  southern  barrier  and  generally  by  shortening  of  the  north¬ 
ern  barrier  (especially  from  1973-1983). 
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The  general  stability  of  a  southern  inlet  location  in  this  system  is 
not  surprising.  Most  of  the  tidal  prism  passes  through  the  deeper  southern¬ 
most  channels  of  the  marsh,  therefore  a  southern  inlet  provides  the  most 
direct  link  to  the  ocean.  Frictional  dominance  in  this  shallow  inlet/ 
estuary  system  (Aubrey  and  Speer,  in  prep.;  Speer  and  Aubrey,  in  prep.) 
makes  this  an  energetically  favorable  location  for  the  inlet.  The  1983 
location  of  the  inlet  requires  that  long  (~2  km)  shallow  channels  carry 
most  of  the  tidal  prism  to  the  south.  A  large  fraction  of  the  total  tidal 
energy  is  dissipated  in  these  channels;  consequently  development  of  an 
energetically  more  favorable  inlet  location  further  south  is  probable  in 
the  near  future  (order  of  a  decade). 

The  southernmost  location  places  constraints  on  barrier  beach  configur¬ 
ation.  In  general,  the  northern  barrier  is  not  strongly  eroded  by  the  ebb 
tidal  flows  when  the  inlet  is  in  a  southerly  location  as  compared  to  a  more 
northerly  one  (reasons  for  this  are  presented  later).  A  short,  slowly- 
growing  southern  spit  can  develop  without  catastrophic  storm  influence. 
However,  large-scale  growth  of  a  southern  barrier  spit  is  dependent  on 
storm  activity  and  breaching  of  the  northern  barrier. 

Three  mechanisms  appear  responsible  for  observed  updrift  tidal  inlet 
migration  at  Nauset  (figure  11):  a)  attachment  of  distal  ebb  tide  delta 
bars  to  the  downdrift  barrier  spit;  b)  storm-induced  breaching  and  sub¬ 
sequent  stabilization  to  form  a  new  inlet;  and  c)  ebb  tide  discharge 
around  the  inlet  channel  bend.  The  last  two  mechanisms  can  result  in 
either  updrift  or  downdrift  migration.  The  first  two  mechanisms  have  been 
observed  widely  at  other  inlets;  the  third  mechanism  is  described  here  for 


the  first  time. 


ACCRETION  OF  EBB  DELTA  BARS  JT  STORM ■  INDUCED  SHIFTS  HI  FLOW  AROUND  A  BEND 
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Figure  11.  Three  modes  of  updrift  inlet  migration  responding  to  different 
combinations  of  waves,  tides  and  storms.  All  three  modes  have 
been  observed  at  Nauset  Inlet,  Cape  Cod,  Massachusetts. 
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a)  Any  model  of  Nauset  Inlet's  migration  must  include  mechanisms  for 
longshore  bypassing  of  sediment  past  the  inlet,  since  the  volume  rate  of 
longshore  sand  transport  is  a  primary  factor  controlling  inlet  stability 
(e.g.,  Brunn  and  Gerritsen,  1959;  Bruun,  1978).  Large  net  longshore  trans¬ 
port  rates,  estimated  to  be  about  250,000  m3/yr  to  the  south  (U.S.  Army 
Corps  of  Engineers,  1969),  occur  at  Nauset  Inlet.  Consequently,  shallow 
overwashes  are  typically  filled  quickly,  and  the  tidal  prism  seems  capable 
of  supporting  only  a  single  stable  inlet. 

A  common  mode  for  sediment  bypassing  of  tidal  inlets  is  through  forma¬ 
tion  and  migration  of  distal  ebb  delta  bars,  which  are  predominantly  wave- 
driven.  For  many  wave-influenced  tidal  inlets  throughout  the  world,  ebb 
delta  bar  migration  is  the  dominant  bypassing  mode  (Bruun,  1978;  Fitzgerald, 
1983;  Nummedal,  1983),  while  bypassing  through  the  inlet  proper  (Galvin, 
1983)  appears  less  important.  Bar  bypassing  results  in  episodic  accretion 
of  sand  on  the  downdrift  barrier,  often  increasing  encroachment  of  the  down- 
drift  barrier  into  the  inlet  throat. 

Bar  bypassing  has  been  documented  several  times  at  Nauset  Inlet 
(figure  12),  increasing  the  length  of  the  southern  (downdrift)  barrier. 

The  existence  of  active  bar  bypassing  at  Nauset  causes  south  spit  to  grow 
to  the  north,  agains t  the  influence  of  predominant  longshore  transport, 
since  the  migrating  bars  weld  to  the  spit  terminus  instead  of  escaping  the 
inlet  influence  (figure  11).  Although  not  shown  on  the  figure,  the 
accreted  sand  remains  on  the  downdrift  spit,  increasing  its  length,  and 
forcing  the  inlet  to  migrate  northwards.  Contrary  to  the  observations  of 
Fitzgerald  (1983),  these  accretionary  episodes  occur  on  time  scales  of 


months,  not  years. 


b)  The  importance  of  storm  activity  to  major  changes  in  inlet/barrier 
beach  configuration  is  illustrated  by  comparison  of  inlet  migration  rates 
with  storm  frequency.  Historical  data  (figures  8,  9  and  10)  show  three 
periods  of  high  storm  activity  since  1933,  preceded  by  a  48  year  period  of 
relative  quiescence.  The  first  period  of  intense  activity  lasted  from  1933 
until  1939.  Unfortunately,  inadequate  chart  and  photo  coverage  prevents 
full  documentation  of  inlet  response  to  this  stormy  period.  The  second 
stormy  period  covered  the  years  1950  to  1962.  Large  scale  inlet  migration, 
together  with  overwash  and  breaching  of  the  barrier  beaches,  occurred  dur¬ 
ing  this  time.  The  north  barrier  breached  in  May  1953  and  January  1956, 
while  the  south  barrier  breached  in  December  1957  and  early  spring  1960. 
Storm-induced  changes  in  barrier  beach  length  of  as  much  as  780  m  have  been 
observed.  A  third  period  of  intense  storm  activity  existed  in  the  early 
1970's.  One  of  the  peak  years,  1972,  coincides  with  a  breach  in  the  north 
spit,  which  initiated  the  current  phase  of  steady  northward  inlet  migration. 

c)  Casual  observation  of  tidal  inlet  and  estuarine  flows  shows  that 
these  channels  often  are  not  straight,  but  rather  have  pronounced  curva¬ 
ture.  This  curvature  has  dramatic  effects  on  flow  through  these  channels, 
affecting  near-bed  shear  stress  distributions  and  resultant  sediment  trans¬ 
port.  Complexity  of  channel  geometry  ranges  from  long,  straight  channels 
with  occasional  bends,  to  nearly  continuous,  sinuous  geometry  reminiscent 
of  river  channel  meanders.  Channel  curvature  within  an  inlet  mouth  pro¬ 
vides  a  mechanism  for  inlet  migration,  as  discussed  below.  An  extensive 
literature  discusses  the  effects  of  channel  curvature  on  flow  structure 
and  sediment  transport,  largely  resulting  from  an  interest  in  river  channel 
meanders  and  open  channel  flow.  Much  of  the  work  to  date  has  been  done  by 
engineers  interested  in  channel  scour  and  deposition  (e.g.,  Nouh  and  Town¬ 
send,  1979),  or  by  geologists  studying  riverine  processes  (e.g.,  Dietrich 
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Figure  13.  Nauset  Inlet  in  summer  of  1980  shows  the  strong  channel  curvature  responsible  for 
'river-bend'  flow  structure.  Primary  tidal  channels  in  the  estuary  are  to  the 
south  (left) ,  causing  this  curved  flow. 
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et  al.,  1982).  Recently,  more  complete  numerical  models  have  quantified 
aspects  of  channel  bend  flow  which  had  been  described  qualitatively  by 
previous  researchers  (e.g..  Smith  and  McLean,  1983). 

Since  the  inlet  mouth  at  Nauset  has  considerable  curvature  (figure  13), 
sediment  transport  patterns  are  modified  by  resulting  bed  shear  stress  gra¬ 
dients.  To  examine  the  magnitude  of  flow  curvature  at  Nauset  Inlet  and 
evaluate  its  effect  on  sediment  transport,  simple  theory  was  developed  and 
tested  by  field  observations  (Aubrey  and  Speer,  1983).  Since  flow  curva¬ 
ture  induces  gradients  in  water  level  through  a  channel  bend,  and  water 
level  is  relatively  easy  to  monitor  compared  to  velocity  or  shear-stress 
distribution,  modeling  and  measurement  efforts  concentrated  on  surface  gra¬ 
dients  across  the  inlet  mouth. 

For  simple  geometry  and  quasi-steady  channel  flows,  the  along-channel  (n) 
momentum  equation  reduces  to  a  simple  balance  between  the  bottom  stress 
(Tb)  and  the  sea  surface  gradient  On/8n): 

xb  =  -pgh  8o/3n  (1) 

where  h  =  mean  water  depth,  p  is  density  of  water,  and  g  is  gravitational 
acceleration.  If  the  downstream  set-down  is  measured,  an  estimate  of  total 
bottom  friction  is  obtained.  The  lowest  order  equation  for  cross-stream 
flow  reduces  to: 

u2  =  8_  (n)  (2) 

gR  8r 

for  the  case  of  small  channel  half-width  (b)  compared  to  radius  of  curva¬ 
ture  (R).  U  is  the  depth-averaged  velocity.  For  a  sinusoidal  channel,  the 
cross-stream  gradient  reduces  to: 

/2b  \ 

(nb-n-b)  =(gR<Ju2 


(3) 


Field  experiments  covering  five  days  (described  in  detail  by  Aubrey 
and  Speer,  1983)  examined  the  magnitude  of  sea  surface  gradients,  and 
checked  their  consistency  with  channel  bend  theory.  Maximum  instantaneous 
down-channel  gradient  over  the  five-day  period  during  ebb  tide  was  0.0007 
(thirty  centimeters  over  the  channel  length  separating  the  two  sensors), 
while  maximum  flood  tide  gradients  were  0.0005  (twenty  centimeters  over 
the  same  separation).  Average  maximum  gradients  over  the  five-day  period 
are  the  same  for  flood  and  ebb  tide  (18.5  cm  over  the  instrument  separa¬ 
tion).  From  these  gradients,  a  maximum  shear  stress  of  approximately  140 
dynes/cm7  was  calculated,  where  these  shear  stress  estimates  include  not 
only  near-bed  friction  but  also  form  drag,  wave/current  interaction,  and 
sediment  transport  effects  (e.g..  Grant  and  Madsen,  1982).  Bed  shear 
stress  is  not  easily  separated  from  this  total  shear  stress.  Maximum  cross¬ 
channel  set-up  of  approximately  5  cm  was  calculated  from  the  down-channel 
gradients,  using  equation  (5),  which  represents  a  slope  of  0.0009,  similar 
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to  the  down-channel  gradient.  Field  measurements  show  a  cross-channel  set¬ 
up  of  complicated  structure,  with  a  magnitude  of  f ive-to-eight  cm,  most  of 
which  is  consistent  with  the  simple  model  presented  above.  The  remainder 
of  the  observed  set-up  is  of  unknown  hydrodynamic  origin  and  not  explained 
by  our  simple  theory. 

Theory  and  observation  are  consistent  with  the  analogy  between  Nauset 
Inlet  flow  structure  and  river-bend  flows,  although  neither  the  measure¬ 
ments  nor  theory  allow  for  in-depth  comparison  of  the  flow  fields.  Result¬ 
ing  inlet  morphology  is  also  consistent  with  the  hypothesis  that  channel 
bends  are  responsible  for  inlet  migration.  Nauset  Inlet  has  both  a  steep 
outer  (northern)  channel  bank,  and  an  accreting  point  bar  on  the  inner 
(southern)  bank,  similar  to  sedimentation  patterns  in  river  bends.  Observ¬ 
ation  of  Nauset  Inlet  since  1972  shows  a  migration  pattern  consistent  with 
the  river  bend  analogy,  with  the  south  spit  elongating  and  the  north  spit 
shortening  due  to  erosion.  Although  analogies  between  inlet  bends  and 
river  bends  based  on  simple  theory  and  observations  of  flow  patterns  and 
sedimentation  are  incomplete,  they  indicate  the  potential  importance  of 
curvature  in  inlet  geometry  on  migration  history  of  tidal  inlets. 

A  conceptual  model  of  inlet  migration  (figure  11)  based  on  the  three 
mechanisms  described  above  explains  the  unusual  northward  movement  of 
Nauset  Inlet,  opposite  the  longshore  transport  direction.  When  the  inlet 
is  at  its  southernmost  location,  the  estuary  channels  empty  directly  into 
the  ocean.  There  is  no  curvature  to  the  flow  (which  could  result  in  com¬ 
plex  flow  non-uniformity),  and  the  north  barrier  is  not  preferentially 
eroded.  Tidal  flows  are  strong  enough,  however,  to  prevent  material  trans- 


ported  past  the  north  spit  from  filling  in  the  inlet  channel.  Bar  bypass¬ 
ing  of  littoral  drift  leads  to  accretion  along  the  downdrift  spit.  If 
storms  do  not  halt  this  growth,  the  south  spit  can  extend  to  the  north  and 
eventually  overlap  the  north  spit,  as  occurred  in  the  late  1960 's  and  early 
1970's.  The  base  of  the  inlet  channel  retains  its  southerly  location,  and 
the  channel  simply  elongates  to  the  north.  In  this  configuration,  the 
north  barrier  remains  essentially  unchanged.  This  particular  pattern  has 
been  observed  once  in  the  last  30  years. 

A  different  mechanism  was  responsible  for  the  barrier  overlap  pattern 
observed  in  the  1950's.  In  this  instance,  the  northern  barrier  was 
breached  during  a  storm.  The  barrier  remnant  south  of  the  breach  attached 
to  Nauset  Heights  to  form  a  relatively  long  southern  spit.  The  base  of 
the  inlet  channel  retained  its  southern  location  and  extended  through  the 
breach.  Subsequently,  the  northern  barrier  lengthened  by  attachment  to  a 
marshy  island  in  the  bay.  The  south  spit  further  elongated  through  bar 
bypassing,  resulting  in  a  pattern  of  overlapping  spits.  As  in  the  1965- 
1972  pattern,  the  inlet  channeL  lengthened  as  the  south  spit  grew  to  the 
north.  The  north  spit  remained  relatively  stable  after  attachment  to  the 
bay  island.  Neither  of  these  patterns  is  representative  of  the  present 
migration  phase.  The  past  ten  years  of  movement  have  been  characterized 
by  actual  northward  migration  of  the  inlet  channel  and  shortening  of  the 
north  spit,  with  no  barrier  overlap. 

A  storm  breaching  the  north  spit  in  1972  caused  this  new  migration 
pattern  to  develop.  The  main  inlet  channel  stabilized  further  north  than 
in  previous  storm  breaches.  As  a  result,  the  dominant  ebb  tidal  flow  was 
constrained  to  flow  to  the  north,  and  then  east  through  the  inlet  channel, 
setting  up  a  channel  bend  flow  pattern  with  erosion  on  the  outer  part  of 
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the  bend  (north),  and  accretion  on  the  inside  of  the  bend  (south).  The 
north  wall  of  the  inlet  channel  is  presently  eroding  while  a  large  sand 
deposit  is  forming  on  the  south  bank  of  the  channel.  This  "flow  around  a 
bend"  has  existed  for  approximately  ten  years,  leading  to  nearly  1  km  of 
northward  inlet  movement.  This  migration  will  probably  continue  until 
either  the  inlet  encounters  an  erosion-resistant  substrate  or  a  major  storm 
changes  the  inlet  location.  In  the  case  of  Nauset,  the  former  is  unlikely 
(see  Aubrey  et  al.,  1982)  since  inlet  tidal  flows  are  currently  eroding 
peat  deposits  underneath  the  sandy  barrier  spits.  Nothing  more  erosion- 
resistant  is  likely  to  be  encountered.  Storm-induced  inlet  relocation  is 
a  strong  possibility.  The  long,  frictionally  dominated  channels  presently 
carrying  the  tidal  prism  would  probably  be  abandoned  if  a  more  southerly 
breach  were  created  by  a  storm.  In  that  case,  the  large  long-shore  trans¬ 
port  could  quickly  close  off  the  present  inlet,  which  is  nearly  clogged  by 
the  extensive  ebb-tide  delta.  A  more  northerly  breach  created  by  storm 
overwash  is  not  likely  to  persist  (as  the  February,  1978  blizzard  demon¬ 
strated)  because  the  northern  depths  and  tidal  prism  are  to'o  small.  An 
additional  factor  increasing  the  likelihood  of  breaching  near  Nauset  Harbor 
is  the  narrow  width  of  the  barrier  at  this  point.  This  narrowing  is  caused 
by  erosion  on  the  bay  side  of  the  barrier  during  ebb  tide,  as  the  easterly¬ 
flowing  tide  is  redirected  northwards  towards  the  present  inlet  (resulting 
in  another  complex,  channel-bend  flow  pattern). 

The  long-term  fate  of  this  estuary  is  affected  by  two  dominant  trends: 
inlet  migration  (which  contributes  sediment  to  the  estuary  via  flood  tide 
delta  growth)  and  westward  spit  migration.  Both  of  these  factors  reduce 
the  tidal  prism,  and  consequently  reduce  the  equilibrium  cross-sectional 
area  of  the  inlet.  The  first  factor  has  been  discussed  in  detail.  The 
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second  factor,  net  onshore  migration  of  the  Nauset  barrier  beach,  is  appar¬ 
ent  in  spite  of  large,  higher  frequency  fluctuations  (Speer  et  al.,  1982). 
The  steady  shoreward  migration  is  a  result  of  sea  level  rise  combined  with 
overwash  and  inlet  processes  (barrier  roll-over).  Higher  frequency  oscil¬ 
lations  superimposed  on  this  steady  retreat  result  from  inlet  migration 
episodes,  seasonal  beach  changes,  bar  bypassing  events  and  large,  near¬ 
shore  bedform  generation  (Aubrey,  1980).  The  effect  of  the  onshore  migra¬ 
tion  is  a  reduction  in  tidal  prism  (specifically  by  reduction  of  the  area 
of  the  back  bay).  Tidal  prism  is  also  reduced  by  deposition  of  sand  as  a 
flood  tide  delta,  an  important  factor  since  1972,  as  the  inlet  has  steadily 
migrated  northwards  approximately  1  km.  Vestiges  of  the  former  flood  tidal 
deltas  are  visible  on  recent  aerial  photographs  of  the  area.  As  a  result 
of  overwash  and  bay  infilling,  the  stable  inlet  configuration  will  become 
narrower  and  shallower  with  time  (reduced  equilibrium  cross-sectional 
area).  This  filling  trend  currently  exceeds  back-barrier  deepening  attri¬ 
butable  to  sea-level  rise,  but  anticipated  increased  rates  of  sea-level 
rise  (Aubrey  and  Emery,  1983)  may  reverse  this  trend. 

SUMMARY 

Three  distinct  patterns  of  natural  inlet  migration  have  been  identified 
from  historical  data  (figure  11),  and  their  underlying  causes  hypothesized. 
These  mechanisms  explain  the  rare  case  where  an  inlet  migrates  in  a  direc¬ 
tion  opposite  the  dominant  longshore  sand  transport,  such  as  at  Nauset 
Inlet.  Large  variability  in  barrier  spit  length  across  a  baymouth  can  also 
be  a  reflection  of  these  mechanisms.  This  rapid,  spatially  variable,  inlet 
migration  contributes  to  infilling  of  some  inlet/estuary  complexes  on  a 
geological  time  scale,  as  the  continually  enlarging  flood  tide  delta 
evolves  at  each  inlet  location.  The  result  is  an  accelerated  shrinkage  of 
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some  estuaries,  with  consequent  reduction  in  inlet  channel  depth  and  width 
(the  decreased  channel  area  responding  to  a  reduced  tidal  prism).  Whether 
or  not  this  flood  tide  delta  growth  significantly  alters  the  fate  of  the 
estuary  depends  on  the  hydraulic  characteristics  of  the  inlet  and  estuary 
(flood  tide  delta  growth  is  a  function  of  flood/ebb  flow  dominance),  as 
well  as  long-term  trends  in  sea-level  rise. 

The  three  distinct  patterns  of  inlet  migration  are  (figure  11): 

1)  Growth  of  the  downdrift  spit  by  addition  of  sediment  from  ebb  tide 
delta  distal  bars:  some  of  these  distal  bars  weld  onto  the  downdrift  spit 
without  escaping  the  inlet  environment.  The  time  scale  of  these  growth 
episodes  is  months,  with  an  associated  spit  growth  on  a  scale  of  100  m. 
Resultant  spit  change  are  relatively  small  compared  to  the  other  two  growth 
mechanisms.  This  mechanism  can  cause  only  updrift  inlet  migration. 

2)  Storm-induced  shifts  in  inlet  position  associated  with  superele¬ 
vated  water  levels:  these  changes  are  rare  but  significant,  with  time 
scales  of  tens  of  years  and  spatial  scales  of  hundreds  of  meters.  Storm 
breaches  will  remain  stable  and  replace  previous  inlets  if  they  are 
hydraulically  more  efficient  than  alternative  breaches.  These  major  inlet 
relocations  have  played  an  important  role  at  Nauset  Inlet,  by  shifting  the 
inlet  position  to  the  north  (against  the  sense  of  net  littoral  drift)  and 
allowing  the  flow  characteristics  to  set  up  a  stable,  steady  northward 
inlet  migration  independent  of  storm  influences. 

Storm  effects  in  the  future  are  expected  to  influence  the  Nauset  barri¬ 
ers  significantly,  and  shift  the  inlet  to  the  south.  Since  the  southern¬ 
most  limit  of  the  estuary/inlet  system  has  historically  been  the  preferred 
position  (because  it  is  the  most  efficient  location  for  tidal  exchanges 
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between  the  ocean  and  bay),  a  breach  at  this  narrow  part  of  the  barrier 
will  likely  become  the  preferred  inlet  position.  At  present,  a  stable 
dune-line  is  inhibiting  storm  overwash  and  breaching  at  this  location. 

3)  Steady  northward  migration  characterized  by  flow  around  a  bend 
(erosion  on  outside  of  bend,  accretion  on  inside  of  bend)  during  ebb  tides 
this  migration  has  occurred  since  1972  when  a  storm  breach  rapidly  shifted 
the  inlet  location,  setting  up  a  long,  confined  southern  barrier-parallel 
channel  through  which  most  tidal  exchange  takes  place.  Ebb  flow  through 
this  barrier-parallel  channel  must  make  a  sharp  bed  through  the  inlet  to 
exit  into  the  ocean.  This  bend  creates  a  distinctive  three-dimensional 
flow  pattern  similar  to  river  bend  flows,  eroding  the  north  spit  and 
accreting  to  the  south.  The  result  is  a  steady  northward  migration  which 
will  cease  when  a  storm  opens  a  breach  further  south  of  the  present  inlet; 
this  new  breach  will  likely  become  the  preferred  inlet  position.  Flow  cur 
vature  resulting  from  complex  inlet/tide  channel  geometry  may  be  respons¬ 
ible  for  both  updrift  and  downdrift  inlet  migration  at  other  locations, 
playing  an  important  role  in  barrier  beach  evolution. 
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CHAPTER  III 

A  STUDY  OF  NON-LINEAR  TIDAL  PROPAGATION 
IN  SHALLOW  INLET/ESTUARINE  SYSTEMS 
PART  I :  OBSERVATIONS 
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ABSTRACT 


The  offshore  tide  becomes  strongly  distorted  as  it  propagates  into 
shallow  estuarine  systems.  Observations  of  sea  surface  elevation  and  hori¬ 
zontal  currents  over  periods  ranging  ,  rom  three  days  to  one  year,  at  nine 
stations  within  Nauset  inlet/estuary,  document  the  non-linear  interaction 
of  the  offshore  equilibrium  tidal  constituents.  Despite  strong  frictional 
attenuation  within  the  estuary,  the  overtides  and  compound  tides  of  M2, 

S2  and  N2,  in  particular,  reach  significant  amplitude,  resulting  in  strong 
tidal  distortion.  High  frequency  forced  constituents  in  sea  surface  are 
phase-locked,  consistently  leading  the  forcing  tides  by  60°-70°,  resulting 
in  a  persistent  distortion  with  falling  tide  longer  than  rising  tide. 

Forced  constituents  in  currents  are  more  nearly  in  phase  with  equilibrium 
constituents,  producing  flood  currents  which  are  shorter  but  more  intense 
than  ebb  currents.  A  compound  fortnightly  tide,  MSr,  modulates  the  mean 
water  level  such  that  lowest  tides  occur  during  Neap  phase  instead  of 
Spring  phase.  This  fortnightly  tide  can  be  contaminated  by  storm  surge, 
changing  the  phase  characteristics  of  this  constituent.  Implications  of 
the  overtides,  compound  tides,  and  lower  frequency  tides  on  near-bed  and 
suspended  material  transport  are  profound. 
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I.  INTRODUCTION 

The  astronomical  tide  is  strongly  distorted  during  its  propagation  from 
offshore  into  shallow  inlet/estuarine  systems  common  to  the  U.S.  East  and 
Gulf  Coasts.  Time  asymmetries  develop  in  the  rise  and  fall  of  the  surface 
tide  with  resulting  time  and  amplitude  asymmetries  in  the  velocity  field 
(e.g.,  Boon  and  Byrne,  1981).  This  distortion  can  be  represented  as  the 
non-linear  growth  of  harmonics  of  the  principal  ocean  astronomical  consti¬ 
tuents  (e.g.,  Dronkers,  1964;  Pingree  and  Griffiths,  1979).  Harmonic 
growth  is  a  result  of  finite  amplitude  effects  entering  through  friction, 
non-linear  advection,  and  interactions  with  channel  geometry  as  the  tide 
oscillates  within  the  estuary.  In  the  inlet/estuarine  systems  of  interest 
to  this  study,  the  water  column  is  well  mixed  throughout  most  of  the  tidal 
cycle  and  freshwater  inflow  forms  a  negligible  part  of  the  tidal  prism. 

The  astronomical  tide  spectrum  is  composed  of  a  large  number  of  consti¬ 
tuents  whose  mutual  non-linear  interactions  represent  a  complex  physical 
problem  (e.g.,  Munk  and  Cartwright,  1966;  Gallagher  and  Munk,  1971). 

Primary  frequencies  of  interest  are  integral  linear  combinations  of  six 
basic  components  related  to  celestial  mechanics  of  the  Earth-Moon-Sun 
system: 

fT1  =  1  day  =  period  of  Earth's  rotation  relative  to  the  Sun 

fi1  =1  lunar  month  =  period  of  Moon's  orbital  motion 

f 3 1  =1  year  =  period  of  Sun's  orbital  motion 

f 4  -  8.85  years  =  period  of  lunar  perigee 

f 5  -  18.61  years  =  period  of  regression  of  lunar  nodes 

ft;1  -  20,900  years  =  period  of  solar  perigee 

Interactions  between  these  basic  frequencies  result  in  other  energetic 

tidal  frequencies,  such  as  the  inverse  lunar  day  (fj1),  where: 


f4  =  f! 


f2  +  f3  =  0.966  fi 
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Tidal  frequencies  are  divided  into  species  which  are  separated  by  one 
cycle/lunar  day,  into  groups  which  are  separated  by  one  cycle/month,  and 
constituents  which  are  separated  by  one  cycle/year.  Here  we  are  concerned 
with  various  predominant  species  and  associated  groups  (table  1). 

Compound  constituents,  which  are  linear  combinations  of  basic  frequen¬ 
cies,  can  be  generated  through  non-linear  celestial  and  fluid  mechanics. 

For  example,  MSf  consists  of  both  a  weak  astronomical  term  and  a  potenti¬ 
ally  larger  hydrodynamic  term  arising  from  M2  -  S2  interactions,  while 
MS4  is  a  constituent  arising  solely  from  non-linear  hydrodynamic  interaction 
of  M2  and  S2.  Similarly  M4 ,  Ms,  and  Mg  have  no  equilibrium  tidal  argument, 
but  reflect  non-linear  generation  in  oceanic  basins  of  various  scales. 
Magnitude  and  phase  of  each  observed  constituent  provides  insight  into 
hydrodynamic  processes. 

The  present  study,  along  with  its  companion  study  (Speer  and  Aubrey, 
1984),  focuses  on  the  distortion  of  the  offshore  tide  as  it  propagates 
through  shallow  estuarine  systems  connected  to  the  ocean  by  a  narrow  tidal 
inlet.  Frictional  decay  of  the  offshore  tide  and  harmonic  growth  of  forced 
constituents  form  the  basis  for  examination  of  non-linear  processes.  Diag¬ 
nostic  numerical  modeling  examines  the  most  important  and  representative 
interactions  observed  in  the  field.  Robinson,  Warren,  and  Longbottom 
(1983)  reported  a  similar,  but  less  comprehensive,  study  conducted  in  the 
Fleet,  on  the  south  coast  of  England. 

Along  the  US  northeast  coast,  the  principal  tidal  constituents  of 
interest  are  the  M2  tide  and  its  harmonics.  The  magnitude  of  M2  results 
in  it  dominating  non-linear  processes  during  tidal  propagation  through  the 
estuary,  although  compound  tides  are  also  generated.  In  particular,  M2  and 
its  first  harmonic,  M4 ,  can  be  used  to  illustrate  the  dominant  features  of 
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TABLE  1 


TIDAL  SPECIES  AND  CONSTITUENTS  OF  INTEREST 


Species 

Constituent 

Fortnightly 

MSf 

1.4 

solar  days 

Diurnal 

Oi 

25.8 

solar  hours 

Ki 

23.93 

solar 

hours 

Semi-Diurnal 

n2 

12.66 

solar 

hours 

m2 

12.42 

solar 

hours 

S2 

12.00 

solar 

hours 

Ter-Diurnal 

mk3 

8.18 

solar 

hours 

Quarter-Diurnal 

mn4 

6.27 

solar 

hours 

m4 

6.21 

solar 

hours 

ms4 

6.10 

solar 

hours 

Sixth-Diurnal 

M6 

4.14 

solar 

hours 

Ss 

4.00 

solar 

hours 
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tidal  asymmetries  in  these  systems  (fig.  1).  Two  (arbitrary)  phase  relation¬ 
ships  between  M2  and  M4  constituents  (sea  surface  or  velocity),  can  be 
defined: 

Am 2  =  ai  cos  (tot  -  0i  ) 

Am 4  =  a2  cos  (2ut  -  02) 

For  the  case  of  sea  surface,  when  M*  leads  M2  by  90°,  falling  tide  exceeds 
rising  tide  in  duration  (fig.  1A).  The  phase  relationship  in  this  situa¬ 
tion  is: 

20i  -  02  =  90° 

If  we  let  the  constituents  refer  to  tidal  velocity,  with  M2  and  M4  in 
phase,  the  case  of  shorter,  enhanced  flood  currents  is  demonstrated.  This 
phase  relationship  can  be  expressed  as: 

20!  =  02 

The  cases  of  symmetrical  and  enhanced  ebb  currents  can  be  shown  in  a 
similar  manner. 

The  development  of  tidal  asymmetries  can  have  important  effects  on  both 
the  geological  evolution  of  shallow  estuaries  as  well  as  on  the  short  term 
navigability  of  estuarine  channels.  An  estuary  characterized  by  shorter, 
more  intense  flood  then  ebb  currents  (flood  dominant)  may  be  unable  to 
flush  entering  sediments  effectively.  Conversely,  an  estuary  with  stronger 
ebb  than  flood  currents  (ebb  dominant)  may  represent  a  more  stable  configu¬ 
ration.  The  magnitude  of  the  velocity  asymmetry  depends  on  the  non-linear¬ 
ity  of  the  tide.  Both  flood-  and  ebb-dominated  inlet/  estuarine  systems 
are  found  along  the  U.S.  east  coast  (e.g.,  Nauset  Inlet,  flood  dominant; 
Wachapreague  Inlet,  ebb  dominant).  In  addition  to  the  different  types  of 
asymmetry,  these  estuaries  also  exhibit  varying  degrees  of  non-linear 
response  to  tidal  forcing  as  measured  by  the  ratio  of  amplitudes  of  M4  to 


TIDAL  VELOCITY  (cm/sec)  TIDAL  ELEVATION  (cm) 
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TIDAL  ASYMMETRIES 

FLOOD  VERSUS  EBB  DOMINANCE 
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Figure  1. 


Examples  of  tidal  asymmetries  for  the  case  of  M4/M2=0.3. 
la)  M4  is  90°  out  of  phase  with  M2 .  lb)  M4  in  phase  and 
180°  out  of  phase  with  M2 . 
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M2  (the  M4/M2  ratio).  To  better  understand  the  problem  of  tidal  propaga¬ 
tion  in  a  shallow  inlet/estuarine  complex,  a  field  program  and  a  numerical 
modeling  study  were  initiated  at  Nauset  Inlet,  Cape  Cod,  MA  (fig.  2).  This 
paper  details  results  of  the  field  experiments  examining  the  tidal  response 
of  this  estuary.  Part  II  (Speer  and  Aubrey,  1984)  focuses  on  diagnostic 
models  of  the  M2,  M4  interactions. 

The  Nauset  Inlet  system  is  a  salt  marsh  intersected  by  three  major 
tidal  channels  and  connected  to  the  ocean  by  a  natural,  unstabilized  inlet 
(fig.  3).  The  offshore  tide  is  predominantly  semi-diurnal  with  a  range  of 
approximately  two  meters.  The  northern  channel  has  extensive  tidal  flats 
and  mean  depths  less  than  one  meter  (fig.  3).  Middle  and  south  channels 
have  regions  of  tidal  flats,  channel  depths  of  2-3  m,  and  terminate  in  a 
deeper  body  of  water.  Town  Cove  (4-6  m).  The  channels  are  well-mixed  over 
most  of  the  tidal  cycle  and  fresh  water  inflow  is  negligible.  Nauset  was 
chosen  for  study  because  in  terms  of  both  channel  geometries  and  tidal 
asymmetry  it  is  representative  of  many  such  systems  on  the  US  east  coast. 
II.  FIELD  METHODS 

Two  field  experiments  were  undertaken  at  Nauset  to  examine  the 
characteristics  of  tidal  propagation  and  to  provide  a  data  set  against 
which  to  test  diagnostic  numerical  models  of  tide  interactions.  A  two- 
week  experiment  in  September  1981  emphasized  measurement  of  currents  near 
the  inlet  channel,  some  month-long  measurements  of  tidal  elevation  through¬ 
out  the  estuary,  and  measurement  of  sea  surface  gradients  through  the  inlet 
proper.  The  second  field  experiment,  extending  from  August-October  1982, 
included  velocity  and  sea  surface  measurements  to  estimate  local  momentum 
balances  and  a  three  month  deployment  of  an  array  of  tide  gages.  Experi¬ 
mental  instrumentation  for  1982  (fig.  3)  consisted  of  Steven's  tide  gauges 
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AGA 


3. 


Location  map  for  Nauset  Inlet  experiment  of  August  through 
October,  1982.  Numbered  locations:  1-Goose  Hummock  (GH; , 
2-Mead’s  Pier  (MP) ,  3-Snow  Point  (SP) ,  4-Middle  Channel 
West  (MCW) ,  5-Middle  Channel  East  (MCE),  6-Nauset  Heights 
(NH) ,  7-Ocean  site,  8-Nauset  Bay  (NB) ,  '  -  9-North  Channel 

(NC).  AGA  refers  to  locations  of  four-day  long  deployments 
of  an  electromagnetic  current  meter  array. 
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(model  nos.  71-A  and  7030,  made  available  by  the  NOAA  National  Ocean 
Survey),  a  five-element  array  of  electromagnetic  current  meters  deployed 
on  a  tripod  (deployments  typically  lasting  2-4  days  with  sampling  rates  of 
1  Hz),  and  numerous  2-7  day  deployments  of  pressure  sensors  (Sea  Data  Model 
TDR-la).  The  inlet  was  instrumented  with  four-day  deployments  of  pressure 
sensors  in  1980  and  1981.  These  latter  measurements  yielded  insight  into 
the  interaction  of  tide  and  the  narrow  inlet  channel,  a  region  where  long 
term  tide  gauging  is  impractical.  Further  details  of  these  experiments 
may  be  found  in  Aubrey  and  Speer  (1983). 

Because  of  the  sensitivity  of  estimates  of  harmonic  tidal  constituents 
to  errors  in  time  base,  attention  was  focused  on  annotating  tide  gage 
records  with  precise  time  codes.  Data  from  these  tide  gages  were  reduced 
to  a  standard  format  either  under  contract  to  NOS  at  their  facilities  (with 
either  a  one  minute  or  six  minute  sample  interval),  or  using  a  Calcomp 
model  9000  coordinate  digitizer  to  convert  analog  records  to  digital  form 
at  the  Woods  Hole  Oceanographic  Institution.  In  both  cases,  proper  quality 
control  on  the  time  base  was  provided.  Other  data  types  were  internally 
recording  on  magnetic  tape,  with  later  processing  at  W.H.O.I. 

III.  STATISTICAL  METHODS 

Statistical  analysis  of  the  data  sought  to  extract  the  amplitudes  and 
phases  of  major  tidal  constituents  from  both  sea  surface  elevation  (tide 
gages  and  pressure  sensors)  and  velocity  (current  meter)  records.  Experi¬ 
mental  logistics  and  instrument  peculiarities  produced  data  records  of 
varying  duration  at  different  areas  of  the  estuary,  making  Fourier  harmonic 
analysis  of  limited  use.  Conventional  fourier  analysis  requires  15-  or 
29-day  record  lengths  to  separate  major  harmonic  components  (M2  from  S2 
and  M2  from  N2,  respectively).  Our  records  were  either  longer  than  this 
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(in  which  case  we  wanted  to  utilize  all  the  data  and  not  just  some  subset) 
or  shorter  than  this,  in  which  case  we  still  wanted  to  extract  tidal 
components . 

For  records  of  15-days  or  more  duration,  initial  analysis  utilized 
standard  power  spectral  techniques  to  determine  the  structure  of  the  tide 
at  various  points  within  the  ocean/inlet/estuary  system.  Vector  velocity 
data  was  analysed  via  rotary  component  spectral  analysis  (Gonella,  1972) 
to  extract  both  frequency  content  as  well  as  any  tendency  for  rotary  tidal 
motion  (negligible  in  long,  straight  channels,  but  potentially  appreciable 
near  channel  bends).  Error  analysis  for  this  method  is  well-established, 
with  power  spectral  values  following  a  chi-squared  distribution,  and  phase 
and  coherence  confidence  limits  described  in  Gonella  (1972)  following  the 
*  work  by  Goodman  (1957). 

Quantitative  estimates  of  harmonic  content  were  derived  from  least- 
square  harmonic  analysis,  a  variant  of  the  standard  Fourier  harmonic 
analysis  technique  developed  by  Schureman  (1971).  The  goal  of  harmonic 
analysis  is  to  extract  amplitudes  and  phases  of  the  primary  tidal  harmonic 
constituents  and  their  compound  tides.  Finite  record  lengths  and  contami¬ 
nation  of  tidal  signals  by  noise  introduced  by  storms  and  other  sources, 
combine  to  make  harmonic  analysis  inexact,  rendering  estimation  of  confi¬ 
dence  intervals  for  these  estimates  of  utmost  importance.  Methods  for  per¬ 
forming  harmonic  analysis  include:  Fourier  Harmonic  Analysis  (Schureman, 
1971),  the  response  or  admittance  method  (Munk  and  Cartwright,  1966; 
Cartwright,  Munk  and  Zetler,  1969),  and  least-squares  harmonic  analysis 
(Boon  and  Kiley,  1978).  The  techniques  vary  in  their  mathematical 
approach,  with  consequent  differences  in  error  estimates.  A  thorough 
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review  of  the  three  techniques  along  with  an  intercomparison  is  provided 
in  Aubrey,  Speer,  and  Boon  (in  prep.);  a  brief  discussion  of  the  techniques 
is  provided  here. 

Fourier  Harmonic  Analysis  is  performed  on  integral  multiples  of  15-  or 
29-day  records  of  sea  surface  elevation  or  velocity.  Record  length  is 
determined  by  the  need  to  resolve  adjacent  energetic  tidal  frequency  bands. 
Spectral  estimates  yield  two  degrees  of  freedom  when  the  minimum  record 
length  for  maximum  resolution  is  used.  These  estimates  have  confidence 
limits  given  by  a  chi-square  distribution,  and  suffer  from  the  internal 
aliasing  due  to  finite  record  length  (gate  function  in  the  time  domain, 
sine  function  in  the  frequency  domain),  since  the  tidal  constituents  are 
not  integral  multiples  of  elemental  frequency  bands.  Low  frequency  and 
unresolved  tidal  constituents  must  be  inferred  using  formulae  such  as  those 
presented  in  bchureman  (1971).  Shorter  record  lengths  are  difficult  to 
quantitatively  examine  using  Fourier  Harmonic  Analysis,  because  of  the 
resolution  restrictions  associated  with  this  method.  In  particular,  separ¬ 
ation  of  the  many  constituents  within  semi-diurnal,  quarter-diurnal,  and 
sixth-diurnal  species  from  a  short  record  is  not  possible.  These  restric¬ 
tions  led  us  away  from  the  Fouiier  Harmonic  Analysis  in  favor  of  other 
schemes . 

The  response  or  admittance  method  was  developed  by  Munk  and  Cartwright 
(1966)  as  an  alternative  to  the  Fourier  Harmonic  Analysis  method.  The 
method  is  based  on  spectrally  relating  a  time  series  of  long  duration  and 
hence  of  excellent  accuracy  to  an  arbitrary  length  series  at  the  same  or 
nearby  location.  The  nomenclature  arises  from  the  admittance  function,  or 
transfer  function,  relating  two  quantities  in  a  spectral  sense.  The 
Fourier  transform  of  the  admittance  function  is  the  impulse  response  func- 
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tion,  which  is  convolved  with  the  input  series  in  the  time  domain  to  yield 
an  output  function.  Alternatively,  the  admittance  (or  transfer)  function 
is  multiplied  by  the  input  transform  to  obtain  the  transform  of  the  output 
function.  Ideally,  the  input  function  is  derived  from  a  suitably  long 
record  of  tidal  data  at  a  reference  station,  generated  either  from  equili¬ 
brium  tidal  theory,  or  lacking  this,  from  a  long  time  series  of  tides  at 
that  coastal  station.  In  the  former  case,  the  admittance  function  yields 
a  relationship  between  the  measured  tide  over  the  period  of  interest  and 
the  equilibrium  tidal  arguments.  Infinite  frequency  resolution  is  poss¬ 
ible  if  the  admittance  function  is  assumed  to  be  smoothly  varying  on  the 
premise  that  (barring  resonance)  the  response  of  the  ocean  to  similar  fre¬ 
quencies  is  similar.  For  harmonic  analysis  purposes,  the  tide  record 
examined  is  related  through  this  admittance  procedure  to  a  noise-free  pre¬ 
dicted  tide  for  the  same  period.  The  part  of  the  measured  tide  that  is 
coherent  with  the  predicted  tide  is  expressed  in  terms  of  a  new  admittance 
function,  whose  statistical  uncertainties  have  been  derived  in  Munk  and 
Cartwright  (1966).  Relative  amplitudes  and  phases  for  the  admittance  func¬ 
tion  derived  from  this  analysis  have  simple  confidence  limits  for  the  case 
of  small  relative  noise  levels. 

The  admittance  technique  is  useful  for  extracting  that  part  of  an 
observed  tide  which  is  coherent  with  the  reference  station  series,  with 
the  remaining  incoherent  part  of  the  signal  interpreted  as  noise.  To  be 
useful,  the  reference  station  and  data  station  must  have  a  similar  response 
to  ocean  tides.  An  admittance  function'  relating  tides  in  the  middle  of  the 
ocean  to  tides  within  an  embayment  has  limited  value,  and  does  not  appreci¬ 
ably  improve  a  prediction  generated  by  an  alternate  harmonic  analysis 
method.  The  extra  work  involved  with  the  admittance  procedure  is  not 
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justified  in  this  situation.  For  the  present  study  of  non-linear  tidal 
harmonic  generation  in  shallow  estuaries,  the  admittance  procedure  was  not 
used  for  several  reasons.  First,  there  is  no  long  time  series  of  tides  in 
the  estuarine  system  from  which  to  extract  response  functions  suitable  for 
generating  noise-free  reference  series.  Second,  no  local  station  is  avail¬ 
able  from  an  environment  physically  similar  to  the  estuarine  system  of 
interest.  A  suitable  reference  station  would  require  similar  non-linear 
response  to  discrete  harmonic  input;  such  a  station  is  not  available.  This 
points  out  the  limited  utility  of  the  admittance  technique  in  many  shallow- 
water  coastal  regions.  The  additional  work  required  for  a  complete  admit¬ 
tance  procedure  represents  another  shortcoming  compared  to  alternate  tech¬ 
niques.  As  a  result,  the  procedure  may  not  be  justified  except  for  certain 
scientific  studies  where  additional  precision  is  made  possible  with  an 
appropriate  long-term  reference  station. 

The  least  squares  harmonic  analysis  method  provides  several  attractive 
features  which  resulted  in  its  adoption  in  this  study.  The  least  squares 
method  can  be  applied  to  any  length  time  series,  extracting  tidal  consti¬ 
tuents  without  the  resolution  restrictions  imposed  by  Fourier  Harmonic 
Analysis,  and  reducing  internal  aliasing  by  extracting  the  exact  consti¬ 
tuent  frequency  from  the  series,  not  the  nearest  integral  multiple  of  the 
fundamental  frequency.  The  nearly  infinite  resolution  of  this  technique, 
termed  super-resolution  by  Munk  and  Hasselman  (1964),  is  attributable  to 
the  least-square  requirement  that  to  determine  phase  and  amplitude  of  a 
signal  of  known  frequency,  only  two  data  points  are  needed.  For  n  frequen¬ 
cies,  then,  only  2n  data  points  are  needed.  No  requirements  are  placed  on 
length  of  record  or  sample  interval.  The  paradox  between  least-squares 
analysis  and  Fourier  harmonic  analysis  is  resolved  by  proper  analysis  of 
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noise  levels  in  the  data  series.  In  particular,  noise  (defined  as  energy 
due  to  non-tidal  sources)  at  frequencies  adjacent  to  tidal  harmonic 
frequencies  limits  the  accuracies  of  harmonic  least  squares  analysis,  such 
that  a  longer  record  will  improve  the  estimates.  Munk  and  Hasselman 
(1964)  first  analysed  the  effect  of  noise  on  least-squares  estimates  of 
two  neighboring  spectral  lines.  Filloux  and  Snyder  (1979)  extended  the 
algorithm  for  the  noise  analysis,  deriving  a  form  for  estimating  phase  and 
amplitude  errors  arising  from  least  squares  analysis. 

Following  Filloux  and  Snyder  (1979),  the  sea  surface  elevation,  C»  at 
tidal  and  lower  frequencies  is  separated  into  harmonic,  Ch»  and  residual, 

Cr,  components: 

C  =  C»  +  C« 

Harmonic  components  are  represented  by  real  variables: 

00 

C«  =  l  nC.  In  =  Cl 
n=l 

and  residuals  become: 

2N 

U  =  C  -  l  n  Ci  In 
n=l 

where  primes  indicated  estimated  quantities,  the  nCi  represents  the 
complex  amplitude  corresponding  to  the  n—  harmonic  frequency,  and  In  can 
be  defined  as: 

I  Zn  -  1  =  COs(unt  +  <$>n  ) 

I  2  n  =  sin(unt  +  4>n) 

4>n  is  the  phase  of  the  corresponding  constituent  of  the  equilibrium  tide 
evaluated  at  time  zero.  The  phase  of  nCi  is  relative  to  the  equilibrium  tide 
at  Greenwich,  equivalent  to  G  in  Schuremann  (1971).  A  scalar  product  is  defined 

t  o  +T 

(a, (3)  =  /  a/3  dt 

to 
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where  T  is  the  record  length.  Conditions  for  solution  of  the  least-squares 
problem  are  set  by  differentiating  (, rZ  with  respect  to  the  nCl  's: 

2N 

l  (In, In.)  mCl  =  (In,C)  n=l  ,  2  ,  .  .  .  2N 
m=l 


Inverted,  this  becomes: 

2N 

mCi  =  I  m„SIl  (I„,C) 
n=l 


where  m„Si  =  (In,Im),  and  Si1  is  the  inverse  of  Si. 
be  decomposed  into  three  parts: 


mC  I  —  1 


+  £  msEl 

^=2N+1 


2N 

a  C I  +  I  nmS I  1  ( In ,£r  ) 

n=l 


These  estimates  can 


N 

Where  m&Ei  =  ^  mnSl  naSl . 

n=l 

The  first  term  is  the  true  complex  amplitude,  the  second  represents  errors 
arising  from  tidal  amplitudes  not  included  in  the  analysis,  and  the  third 
term  represents  errors  due  to  residual  fluctuations.  The  second  term  can 
be  reduced  by  including  a  sufficient  number  of  constituents  in  the  analysis; 
however  these  additional  constituents  adversely  affect  the  third  term.  This 
study  includes  twenty-nine  terms  in  the  analysis,  with  only  four  (MSf,  M2,  M4 , 
M6 )  of  primary  interest.  Consequently,  the  second  error  term  is  small;  the 
third  term  is  the  larger  error  source.  Using  a  definition  of 

CD 

Er  (w)  =  1  J  dt  C  (t)cosut 

2ir  Cr 

as  the  two-sided  spectrum  of  £r,  calculated  from  the  covariance: 

C  (t)  =  <Cr  ( x , t )  ^ R ( x , t +x ) > 

Ci 

where  angle  brackets  represent  ensemble  averaging,  a  co-variance  matrix  can  be 


calculated, 
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<nCi  M>  -  <nCi>  <mCl> 

where  the  summation  is  over  adjacent  constituents  within  the  same  tidal 
species.  Using  appropriate  assumptions  of  quasi-stationarity  in  £r  and 
length  of  record  T,  Filloux  and  Snyder  (1979)  show: 

<nCi  mCi>  -  <nCi>  <n.C{> 

—  ITTEr  (tOn  npSl  mpSl 
P 

where  p,m,n  refer  to  constituents  of  the  same  species.  If 

n'l’H  =  _arg  „£h  =  “tan  1  ( 2n- i£  i /znC  i  ) 

|  n  C  H  |  =  mod  n  C  H  =  ( 2  n - l C I  +  2n£{2)1/2' 

and  a  covariance,  V,  is  defined  for  brevity 
V(a,B)  =  <aB>  -  «x>  <B> 
then  the  variances  of  |  „Ch  |  and  n4»H 
are  given  as: 

I  n^H  I  2  V(  |  n  C  H  I  ,  I  n^H  |  )  = 

2n-lCl2  V(2n-lCl*  2n-l£l) 

+  2  2  n  -  l  ^  1  2 n£ I V( 2 n -  1 C I  »  2nO) 

+  2nCl2  V(2nCl»  2nCl) 

|n^H|4  V(n(t>H,n$H)  = 

2p-lCl2  V(2n-lCl»2n-lCl) 

-  2  2  n -  1 C I  ZnCl  V(zn-lCl»  Zn£l) 

+  2nCl2  V(2n£l»2nCi) 

Errors  can  be  presented  as  standard  deviations,  or 
V (  |  n  C  H  |  ,  (  nC  H  |  )  '  2  and  V(n$H»  n$H  )  *  7 2  . 

Although  algebraically  complex,  these  expressions  show  the  direct  dependence 
of  errors  on  the  noise  spectrum,  Er,  while  the  Aw  dependence  shown  by 
Munk  and  Hasselman  (1964)  is  implicit  in  the  matrix  Si1. 
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For  this  study,  twenty-nine  tidal  constituents  were  extracted  in  a 
least-squares  sense,  fourteen  calculated  directly  with  the  remainder  infer¬ 
red  using  formulae  provided  in  Schureman  (1971).  Of  these,  six  species 
were  of  particular  interest  (table  1),  with  all  species  except  diurnal  and 
semi-diurnal  containing  considerable  forced  constituents.  Once  these 
constituents  were  eliminated,  the  residual  signal  (regenerated  from  the 
original  series  less  the  harmonic  components)  was  analysed  with  a  fast 
Fourier  transform,  to  generate  estimates  of  the  noise  spectrum,  Er.  These 
noise  estimates  were  used  to  calculate  errors  in  phase  and  amplitude 
presented  in  tables  2,  3,  and  4. 

IV.  RESULTS  AND  DISCUSSION 

Harmonic  analyses  of  29-day  records  of  sea  surface  elevation  indicate 
that  Nauset  is  a  strongly  non-linear  as  well  as  frictionally  dominated 
inlet/estuarine  system.  The  non-linear  response  to  tidal  forcing  is 
reflected  in  growth  of  high  frequency  Mz  overtides,  compound  tides  and  a 
forced  low  frequency  MSf  constituent  throughout  the  estuary.  The  frictional 
nature  of  the  estuary  is  demonstrated  both  by  large  phase  changes  in  the  tide 
(~76°  in  2  km  for  the  M2  constituent)  and  by  amplitude  decay  of  the  total 
spectrum  (as  much  as  57  percent  of  the  ocean  range  in  2  km). 

Nine  stations  provide  tidal  data  for  examination  of  tidal  distortion  in 
Nauset  estuary  (fig.  3).  Presentation  of  these  results  is  divided  into 
three  sections:  a)  decay  of  the  primary  diurnal  and  semi-diurnal  constitu¬ 
ents,  and  generation  of  the  overtides  of  M2  (table  2);  b)  generation  of  com¬ 
pound  tides  at  frequencies  gi°ater  than  semi-diurnal  (table  3);  and  c)  genera¬ 
tion  of  the  forced  MSf,  or  fortnightly,  compound  tide  (tables  3  and  4). 


TABLE  2:  FUNDAMENTAL  TIDES  AND  OVERTIDES 
(amplitude  in  cm,  Greenwich  phase  in  degrees, 
confidence  limits  of  one  standard  deviation) 
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TABLE  2:  FUNDAMENTAL  TIDES  AND  OVERTIDES  (continued) 
(amplitude  in  cm,  Greenwich  phase  in  degrees, 
confidence  limits  of  one  standard  deviation) 
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TABLE  3:  COMPOUND  TIDES 
(amplitude  in  cm,  Grenwich  phase  in  degrees, 
confidence  limits  of  one  standard  deviation) 

Station  t0  MSf  MK3  MN4  MS4 


a 

0 

a 

0 

a 

0 

a 

0 

29/09/82 

7.8 

214 

0.1 

107 

0.1 

262 

0.5 

290 

(0TG)  7 

30/10/82 

4.3 

59 

0.1 

349 

0.3 

313 

0.2 

331 

23/09/82 

10.4 

254 

1.4 

101 

2.3 

89 

2.7 

162 

+0.7 

+27 

+0.2 

+6 

+0.2 

+5 

(MCE)  5 

06/10/82 

9.0 

297 

1.6 

96 

2.3 

92 

2.8 

157 

+0.6  +23  +0.3  +7  +0.3  +6 


(MCW )  4  18/10/82  7.3  29  2.0  126  2.8  113  2.7  173 

+0.8  +22  +0.3  +6  +0.3  +6 


30/09/82 

12.1 

255 

1.5 

+0.7 

86 

+28 

2.6 

+0.3 

53 

+7 

3.0 

+0.3 

125 

+6 

18/10/82 

7.5 

30 

2.1 

119 

3.1 

106 

2.7 

164 

+0.9  +24  +0.3  +6  +0.3  +7 


25/08/82 

6.0 

33 

1.5 

+0.8 

182 

+32 

2.8 

+0.3 

110 

+6 

2.4 

+0.3 

210 
+  7 

23/09/82 

(SP)  3 

9.5 

254 

1.7 

+0.8 

152 

+27 

2.8 

+0.4 

122 

+7 

3.0 

+0.4 

199 
+  7 

18/10/82 

7.7 

32 

2.3 

+0.8 

150 

+20 

3.0 

+0.4 

138 

±8 

2.9 

+0.4 

190 

±9 

TABLE  3:  COMPOUND  TIDES  (continued) 
(amplitude  in  cm,  Grenwich  phase  in  degrees, 
confidence  limits  of  one  standard  deviation) 


Station  tQ  MSf  MK3  MN4  MS4 


a 

e 

a 

0 

a 

e 

a 

e 

(NC) 

9 

13/09/81 

3.0 

338 

1.4 

178 

3.7 

146 

3.6 

207 

25/08/82 

6.2 

31 

2.2 

237 

3.4 

161 

3.1 

264 

+1.1 

+28 

+0.5 

+8 

+0.5 

+9 

(NB) 

8 

18/10/82 

9.5 

37 

3.2 

200 

3.6 

178 

3.7 

224 

+  1.4 

+25 

+0.8 

+12 

+0.8 

+12 

17/03/83 

8.4 

44 

1.3 

176 

4.5 

121 

3.6 

235 

(MP) 

2 

26/06/83 

6.3 

71 

2.2 

208 

2.8 

192 

0.9 

351 

22/07/82 

6 . 6 

54 

2.7 

218 

3.7 

150 

3.5 

255 

(GH)  1 


17/03/83  8.3  48 


1.5  184 


5.9  114 


5.1  235 
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TABLE  4  MSf  TIDAL  DATA 


STATION/DATE 

25/8/82 

30/9/82 

18/10/82 

(OTG) 

7 

Amplitude  (m) 

0.043***** 

0.078 

%Total  Tide 

0.2 

0.6 

Relative  Phase 

-37° 

177° 

(MCE) 

5 

Amplitude  (m) 

0.104 

0.090 

%  Total  Tide 

2.0 

1.5 

Relative  Phase 

141° 

4° 

(MCW  ) 

4 

Amplitude  (m) 

0.073 

%  Total  Tide 

1.3 

Relative  Phase 

-4° 

(NH) 

6 

Amplitude  (m) 

0.121 

0.075 

%  Total  Tide 

3.0 

1.3 

Relative  Phase 

140° 

-4° 

(SP) 

3 

Amplitude  (m) 

0.060 

0.095 

0.077 

%  Total  Tide 

1.1 

2.2 

1.7 

Relative  Phase 

12° 

141° 

-9° 

(NC) 

9 

Amplitude  (m) 

0.030**** 

%  Total  Tide 

0.3 

Relative  Phase 

60° 

(NB) 

8 

Amplitude  (m) 

0.062 

0.095 

%  Total  Tide 

1.8 

3.8 

Relative  Phase 

7° 

-26° 

(MP) 

2 

Amplitude  (m) 

0.084* 

0.063*** 

%Total  Tide 

1.6 

1.0 

Relative  Phase 

9° 

7° 

(GH) 

1 

Amplitude  (m) 

0.083* 

0.066** 

X  Total  Tide 

1.5 

1.2 

Relative  Phase 

5° 

5° 

* 

Starts 

17/3/83 

** 

Starts 

22/7/82 

*** 

Starts 

26/6/83 

***  * 

Starts 

13/9/81 

***** 

Starts 

30/10/82 

Table  4 

-77- 


Nonlinearity  in  the  tides  can  be  illustrated  by  comparing  the  offshore 
tide  with  the  estuary  tide  (figs.  4a  and  4b,  table  2).  The  equilibrium 
diurnal  and  semi-diurnal  tides  decay  in  their  propagation  through  the 
estuary,  both  through  friction  and  through  non-linear  transfers  of  energy. 
Non-linear,  forced  tides  are  negligible  in  the  ocean  station,  but  are 
significantly  enhanced  within  the  estuary.  In  terms  of  sea  surface  eleva¬ 
tion,  phase  relationships  between  the  equilibrium  and  forced  constituents 
cause  a  definite  asymmetry  in  estuary  tides,  with  rising  tide  shorter  than 
falling  tide  (fig.  5).  Duration  and  magnitude  asymmetries  also  exist  in 
velocity  records,  with  stronger  flood  currents  than  ebb  currents. 

The  dominant  diurnal  tides  are  the  Ki  and  Oi  constituents,  with  roughly 
equivalent  amplitudes  (fig.  4).  M2  is  the  dominant  semi-diurnal  constituent, 
with  an  amplitude  six  times  greater  than  both  the  N2  and  S2  constituents. 
These  relative  amplitudes  caused  us  to  focus  theoretical  modeling  on  over¬ 
tides  of  M2,  rather  than  overtides  of  N2  or  S2,  or  compound  tides  between 
M2  and  either  S2  or  N2  (e.g.,  MN-* ,  MS.*),  although  the  physics  of  inter¬ 
actions  are  identical.  The  principal  overtide  is  M<* ,  with  smaller  forced 
constituents  arising  from  MK3,  MN<* ,  and  MS.* .  The  forced  MSf  (fort¬ 
nightly)  tide  is  the  largest  compound  constituent. 

a)  Decay  of  diurnal  and  semi-diurnal  tides,  and  generation  of  M2  over¬ 
tides:  The  behavior  of  the  diurnal  and  semi-diurnal  tides  as  they  propagate 

through  the  estuary  reflects  dissipation  due  to  friction  and  non-linear 
spectral  energy  transfer.  Friction  reduces  the  amplitude  of  the  tide  and 
creates  phase  lags,  while  other  non-linear  transfers  reduce  tidal  amplitude 
with  characteristic  phase  change  signatures. 

The  frictional  nature  of  the  estuary  is  demonstrated  by  the  significant 
decay  of  all  semi-diurnal  and  diurnal  constituents,  with  associated  large 
phase  lags  (table  2  and  fig.  6).  The  amplitude  decay  of  M2  from  the  ocean 


OCEAN  TIDES 

NAUSET  INLET.  MA 
30/10/82-28/11/82 


<0 


Figure  4a.  Periodogram  of  offshore  tides  (10  m  water  depth)  outside  of 

Nauset  Inlet,  showing  selected  (out  of  35)  tidal  constituents 
calculated  with  least-squares  harmonic  analysis  procedures. 
Overtides  and  compound  tides  are  clearly  negligivle  here. 
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Figure  4b.  Periodogram  of  sea  surface  tides  within  Nauset  Inlet  at  Nauset 

Bay,  showing  the  distinct  growth  of  overtides  and  compound  tides, 
compared  with  the  offshore  tides  of  figure  4a.  Only  selected 

(10  of  35)  constituents  are  depicted.  Both  figures  4a  and  4b  are 
estimates  based  on  29-days  of  data. 


NAUSET  INLET  TIDES 

TIDAL  ASYMMETRY 


Figure  5.  The  effect  of  tidal  distortion  on  sea  surface,  shows  as  a  24-hour 
reconstruction  of  tides  consisting  of  the  M2  constituent  and  its 
first  over tide,  M4,  with  their  proper  phase  relationships.  The 
offshore  tide  is  nearly  a  pure  sinusoid,  whereas  the  Nauset  Bay 
tide  is  distorted  by  over tides,  leading  to  shorter  flood  tides 
and  longer  ebb  tides. 


M2  PHASE  (deg)  AMPLITUDE  (m)  M4/M2  RATIO 

50  SO  ID  140  0.8  12  16  2.0  0.0  0.1  0.2 
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TIDAL 
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Figure  6.  Tidal  distortion  changes  as  a  function  of  distance  into  the 

estuary  for  each  of  the  primary  drainage  channels.  Top  panel 
shows  sea  surface  M4/M2  ratios;  middle  panel  shows  tidal  ampli¬ 
tude  and  bottom  panel,  M2  phase.  Tide  gage  stations  are  indi¬ 
cated  on  the  abscissa. 
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to  the  far  reaches  of  the  estuary  is  approximately  60%  to  the  north  and  50% 
to  the  south.  Tidal  dissipation  occurs  at  different  rates  along  different 
reaches  of  the  estuary  (fig.  6),  and  is  largest  in  shallow  channels  and  in 
channel  reaches  where  tidal  flats  are  developed.  Along  middle  channel, 
the  greatest  decline  in  amplitude  is  in  the  short  (0.7  km)  reach  between 
stations  5  fMCE)  and  b  (MCW).  This  reach  is  characterized  by  large  tidal 
flats  and  extremely  variable  channel  geometry  over  a  tidal  cycle.  Ampli¬ 
tude  decay  down  the  rest  of  the  channel  to  Town  Cove,  a  reach  characterized 
by  less  time  variable  geometry  and  no  large  regions  of  tidal  flats,  is  much 
less.  North  channel,  with  extremely  shallow  depths  and  extensive  flats  has 
the  highest  rate  of  M2  and  total  tide  amplitude  decay.  Phase  changes  in  M2 
(table  2  and  fig.  6)  display  the  same  pattern  as  frictional  decay  with 
greatest  lags  occurring  in  shallow  channels  and  channels  bordered  by  tidal 
flats . 

The  pattern  of  tidal  dissipation  shows  the  estuary  acts  as  a  low  pass 
filter,  with  larger  decay  rates  for  semi-diurnal  than  diurnal  constituents. 
At  Goose  Hummock  (station  1),  the  semi-diurnal  energy  (averaged  over  N2,  Mz 
and  S2)  is  only  27%  of  that  at  the  ocean  gage,  whereas  the  diurnal  energy 
(averaged  over  Ki ,  0i)  is  43%  of  that  at  the  ocean  gage.  At  Nauset  Bay 
(station  8),  the  semi-diurnal  energy  (averaged  over  N2,  M2  and  S2)  is  only 
21%  of  that  at  the  ocean  gage,  whereas  the  diurnal  energy  (averaged  over  Ki, 
Oi )  is  36%  of  that  at  the  ocean  gage.  From  a  phase  standpoint,  the  semi¬ 
diurnal  tides  have  larger  lags  through  the  estuary  than  do  diurnal  tides. 
Because  of  non-linear  growth  of  higher  frequency  compound  tides,  and  lack 
of  suitable  directly  calculated  (not  inferred)  lower  frequency  elementary 
tidal  constituents,  the  complete  filter  characteristics  (cutoff  frequency, 
reject-band  dropoff  rate,  phase  behavior,  etc.)  cannot  be  defined. 
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Non-1 inear  growth  of  the  harmonics  of  M2  follows  a  consistent  pattern 
through  the  estuary.  The  primary  indicato  *  of  non-linearity  used  here  is  the 
M4/M2  ratio,  which  reflects  the  combined  effects  of  spectral  energy  trans¬ 
fer  from  M2  to  M4  and  frictional  dissipation  through  the  estuary.  As 
discussed  earlier,  the  estuary  tends  to  act  as  a  low  pass  filter  in  the 
absence  of  non-linear  effects;  the  nonlinear  effects  selectively  enhance 
specific  high  frequency  (>2cpd)  compound  constituents.  If  nonl inearities 
were  absent,  therefore,  the  ocean  magnitudes  of  M4 ,  M6,  and  Mg  would 
rapidly  decrease  within  the  estuary  from  their  small  offshore  values.  In 
fact,  the  general  pattern  in  all  three  channels  is  a  steady  increase  in  the 
amplitude  of  M4  with  distance  into  the  estuary.  The  amplitude  of  M4  with¬ 
in  the  estuary  can  be  contrasted  with  the  offshore  tide  (table  2)  which  has 
an  M4  amplitude  less  than  1  cm.  The  rate  of  growth  of  M4  differs  in  the 
three  channels.  In  a  distance  of  2  km  through  north  channel,  where  M4 
reaches  11  cm  in  amplitude,  the  M4/M2  ratio  reaches  0.27.  Middle  channel, 
with  extensive  tidal  flats  but  greater  mean  depths  than  north  channel,  has 
a  slower  rate  of  non-linear  growth.  M4  has  an  amplitude  of  7  cm  at  a  dis¬ 
tance  of  1.7  km  from  the  inlet  mouth.  The  slowest  rate  of  non-linear  growth 
is  found  along  south  channel  where  M4  attains  an  amplitude  of  7  cm  in  2.4  km. 
Proceeding  along  middle  channel  into  Town  Cove,  M4  grows  steadily  to  11  cm 
a  distance  of  6  km  from  the  inlet.  At  this  location,  M4  is  the  second 
largest  constituent  behind  M2.  The  higher  harmonics,  M6  and  M8  ,  are  not  an 
important  part  of  the  sea  surface  spectrum,  with  amplitudes  typically  1  cm 
or  less.  They  do  not  dispLay  any  systematic  pattern  of  growth  through  the 
estuary,  a  result  of  both  limited  tide  gauge  accuracy  and  the  large  high- 
frequency  attenuation  characteristics  of  the  estuary. 
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Throughout  the  estuary,  M4  is  phase  locked  to  M2  with  M4  leading  by 
63°  (standard  deviation  of  1.3°).  The  largest  deviation  from  this  phase 
relationship  is  at  Nauset  Heights,  located  in  a  shallow  embayment  off  the 
main  channel  (fig.  3),  where  M4  leads  M2  by  71°.  The  result  of  this  pat¬ 
tern  of  increasing  M4  amplitude  and  M4  phase  lead  of  63°  is  a  distorted 
surface  tide  with  falling  tide  exceeding  rising  tide  in  duration  (fig.  5). 
The  duration  asymmetry  due  tc  che  dl2,  M4  interaction  becomes  more  pro¬ 
nounced  far  from  the  inlet  and  in  shallow  regions  with  extensive  flats.  It 
reaches  extreme  values  in  Town  Cove  and  in  north  channel  where  ebb  averages 
more  than  7.5  hours  in  duration  (Table  5). 

The  non-linear  spectral  transfer  from  M2  to  M4  accounts  for  only  a 
small  fraction  of  the  total  decay  in  M2,  the  larger  portion  attributable  to 
frictional  dissipation.  At  Goose  Hummock  (station  1),  potential  energy  con¬ 
tained  in  M4  represents  only  about  2%  of  the  total  decrease  in  Mz  potential 
energy  from  the  ocean.  Although  the  transfer  to  is  important  for  water 
and  sediment  transport  purposes,  it  is  not  the  major  contributor  to  decay  of 
M2. 

Despite  gradual  changes  in  the  geometry  of  the  estuary  caused  by  north¬ 
ward  movement  of  the  inlet  (Aubrey  et  al.,  1984),  the  monthly  averaged 
pattern  of  tidal  propagation  appears  consistent  over  a  1-2  year  period 
(Table  2).  Data  from  Snow  Point  (station  3)  demonstrate  similar  phase 
relationships,  amplitude  decay  and  non-linear  growth  despite  the  occurrence 
of  a  large  storm  during  October,  1982.  If  one  examines  individual  neap/ 
spring  cycles  at  five  stations  within  a  month,  however,  the  non-linear 
distortion  of  the  tide  varies  considerably.  During  spring  tides,  the 
increased  tidal  amplitude  due  to  the  in-phase  nature  of  S2  and  M2  enhances 
finite  amplitude  effects  on  the  propagating  tide.  This  is  indicated  by  the 


TABLE  5 


TIDAL  DURATION  ASYMMETRY  AT  NAUSET 

STATION  DISTANCE  FROM  INLET  AVERAGE  EBB  AVERAGE  FLOOD 


Ocean 

— 

6 

hrs . 

13 

min. 

6 

hrs . 

13 

min 

Inlet 

— 

6 

hrs . 

30 

min. 

5 

hrs . 

57 

min 

North  Channel 

1.2  km 

7 

hrs . 

24 

min. 

5 

hrs . 

01 

min 

Nauset  Bay 

2.0  km 

7 

hrs . 

43 

min. 

4 

hrs . 

42 

min 

Middle  Channel 

East 

0.99 

6 

hrs . 

46 

min. 

5 

hrs . 

39 

min 

Middle  Channel 

West 

1.7 

7 

hrs. 

00 

min. 

5 

hrs . 

25 

min 

Nauset  Heights 

2.4 

7 

hrs . 

03 

min. 

5 

hrs . 

22 

min 

Snow  Point 

3.0 

7 

hrs . 

12 

min. 

5 

hrs . 

13 

min, 

Goose  Hummock 

6.1 

7 

hrs . 

32 

min. 

4 

hrs . 

53 

min 
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change  in  the  ratio,  a/h,  where  a  =  amplitude  of  semi-diurnal  tide;  h  = 
mean  water  depth.  Typically,  the  amplitude  of  the  semi-diurnal  tide  is  1.5 
times  larger  at  spring  than  neap  tide.  As  a  result  of  enhanced  non-linear 
effects  during  spring  tide,  the  M4  constituent  is  2.5  to  5  times  larger 
than  at  neap  (Table  6),  such  that  at  some  stations  the  time  asymmetry  of 
the  tide  virtually  disappears  during  neap  tides  (fig.  7).  This,  in  turn, 
should  cause  considerable  variability  in  patterns  of  sediment  transport 
over  the  spring/neap  cycle. 

The  strong  non-linear  effects  visible  in  the  surface  tide  are  also 
present  in  tidal  currents  within  the  estuary.  Spectral  energy  transfer 
is,  in  fact,  more  pronounced  in  the  velocity  field  with  friction  directly 
removing  energy  from  the  primary  semi-diurnal  constituent.  A  velocity 
spectrum  from  a  week-long  current  meter  deployment  in  south  channel  (1.5  km 
from  inlet)  illustrates  the  content  of  tidal  velocities  (fig.  8).  This 
harmonic  generation  particularly  of  the  quarter-diurnal  constituent,  again 
results  in  a  distortion  of  the  tide.  The  type  of  distortion  is  dependent 
on  the  relative  phasing  of  semi-  and  quarter-diurnal  constituents.  Unlike 
the  case  of  the  surface  tide,  M2  and  its  first  harmonic  are  more  nearly  in 
phase.  Considerably  shorter  velocity  records  prevent  us  from  defining  the 
phase  as  carefully  as  for  the  surface  tide. 

The  in-phase  nature  of  M2  and  M<?  produces  a  velocity  signal  generally 
characterized  by  more  intense  flood  than  ebb  currents.  Depending  on 
channel  geometry,  this  effect  becomes  pronounced  further  in  the  estuary 
where  the  tide  is  more  non-linear  (fig.  9).  Continuity  arguments  support 
the  observation  that  a  long  ebb  flow  must  be  balanced  by  a  shorter,  more 
intense  flood  flow.  This  simple  picture  is  complicated  in  certain  regions 
of  the  estuary  by  large  changes  in  channel  geometry  which  occur  over  a 
tidal  cycle.  Ebb  velocities  in  South  Channel  exceed  flood  velocities 
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TABLE  6  Spring/Neap  Response 
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Figure  7.  Differences  in  harmonic  growth  and  tidal  distortion  are  shown  for 
Spring  and  Neap  cycles  (4-day  periods  within  each) .  During  neap 
tides,  tidal  distortion  is  much  less,  as  flood  and  ebb  duration  are 
approximately  equal.  In  Spring  tides,  flood  tide  is  markedly 
shorter  than  ebb. 
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Figure  8.  Periodogram  (using  Fourier  analysis)  of  tidal  velocities  in  South 
Channel  for  a  week-long  period.  fj  is  the  M2  tide,  with  over tides 
(and  compound  tides  which  are  not  resolved  with  this  short  record) 
shown  at  higher  frequencies.  Velocity  is  clearly  more  non-linear 
than  sea  surface. 
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Figure  9.  Three-day  records  of  velocity  (top  panel-taken  2.9  m  above 
bottom  in  3.5  m  water  depth)  and  sea  surface  at  Mead's 
Pier  (location  2,  figure  3). 
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despite  the  time  asymmetry  favoring  ebb.  This  pattern  occurs  because  peak 
ebb  flow  is  confined  to  a  much  smaller  channel  at  this  location. 

b)  Generation  of  compound  tides  at  frequencies  greater  than  semi¬ 
diurnal  other  than  pure  harmonics  of  M2:  Tides  resulting  from  interactions 
of  various  diurnal  and  semi-diurnal  constituents  were  examined.  Although 
results  show  pure  harmonics  of  M2  to  be  the  strongest  overtides  generated 
within  Nauset  Inlet,  other  forced  tides  are  energetic.  The  compound  tides 
discussed  here  include  ter-diurnal  (Mks)  and  quarter-diurnal  (MN*  and  MS4). 
Additional  forced  tides  at  these  and  other  frequencies  exist,  but  were 
smaller  than  those  listed. 

The  rationale  for  examining  only  specific  forced  tides  and  not  the  full 
spectrum  of  possible  constituents  arises  from  both  observation  and  theory. 
Fourier  analysis  of  residual  tide  data  indicates  the  compound  tides  select¬ 
ed  for  analysis  by  the  least  squares  procedure  contain  most  of  the  energy 
in  any  particular  species.  More  constituents  can  be  selected  because  of 
the  "super  resolution"  characteristics  of  least-squares  analysis.  However, 
the  low  residual  energy  remaining  after  removal  of  the  major  compound 
constituents  indicates  additional  ones  would  not  have  statistical  signifi¬ 
cance  at  Nauset.  The  tides  selected  in  this  study  represent  first  order 
non-linear  interactions  between  equilibrium  constituents  with  largest  rela¬ 
tive  energy  levels.  In  fact,  as  Munk  and  Gallagher  (1971)  demonstrate, 
each  of  the  forced  frequency  bands  includes  a  large  number  of  constituents 
involving  2nd  and  3rd  order  interactions  between  forced  and  equilibrium 
constituents.  Non-linear  processes  produce  strong  frequency  splitting  in 
the  non-as tronomica 1  or  forced  bands.  The  strength  of  these  constituents 
in  the  resulting  spectrum  depends  a  great  deal  on  the  relative  energy  level 
of  the  constituents  from  which  they  were  formed.  Energy  levels  are  low  for 
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all  but  five  astronomical  frequencies  (M2,  S2,  N2,  Oi,  Ki)  and  drop  off 
rapidly  past  the  quarter-diurnal  band.  This  fact  along  with  the  dissipative 
tendency  of  Nauset  to  attenuate  higher  frequency  components  of  the  spectrum 
explains  why  most  of  the  residual  energy  is  found  in  relatively  few  compound 
tides . 

The  ter-diurnal  tide  (MK3)  arises  from  interaction  of  M2  and  Ki  consti¬ 
tuents.  Other  possible  ter-diurnal  constituents  not  considered  here  include 
2MKi  (arising  from  2M2-K;),  SK3  (arising  from  S2  and  Ki),  and  SC?,  (arising 
from  S2  and  Oi).  Offshore,  MK3  is  negligible  (table  2),  but  once  inside  the 
inlet  it  increases  progressively  from  a  value  of  1.5  cm  near  the  inlet,  to 
a  maximum  of  3  cm  at  the  far  reaches  of  the  estuary.  It  is  roughly  a  factor 
of  three  less  in  amplitude  than  the  M4  overtide.  The  phase  of  the  forced 
MK3  tide  uniformly  leads  the  phase  of  the  sum  of  the  M2  and  Ki  tides  by  101° 
(+  15°  standard  deviation),  so  it  too  is  phase  locked  to  the  forcing  func¬ 
tion,  although  not  as  narrowly  as  the  over tide.  Another  ter-diurnal 

constituent  of  equal  magnitude  is  MO3,  resulting  from  the  (M2+0i)  interac¬ 
tion.  This  was  not  examined  in  detail  because  its  response  was  similar  to 
MK3.  The  MK3  constituent  does  not  produce  time  asymmetry  in  the  surface 
tide,  but  does  have  a  minor  effect  on  the  diurnal  inequality. 

The  quarter-diurnal  compound  tides  include  the  major  constituents  MN4 
and  MS4,.  which  are  derived  from  interactions  of  (M2+N2)  and  (M2+S2),  respec¬ 
tively.  Other  quarter-diurnal  tides  exist,  including  MK4  (M2+K2),  and  S4 
(2S2),  these  have  not  been  considered  because  of  arguments  similar  to  those 
of  Gallagher  and  Munk  (1971).  MN4  and  MS4  are  approximately  of  equal  magni¬ 
tude,  increasing  uniformly  from  a  value  of  about  2.3  cm  near  the  inlet 
proper,  to  about  3.5  cm  near  the  extremes  of  the  estuary.  The  amplitudes  of 
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both  are  negligible  offshore.  Although  non-linear  growth  is  slower  than  for 
M4 ,  reaching  only  about  one-tenth  its  potential  energy  in  the  far  reaches  of 
the  estuary,  non-linear  transfer  into  compound  tides  exceeds  frictional 
decay.  Phasing  of  these  two  constituents  follows  that  of  the  other  compound 
tides  and  overtides.  MN4  leads  the  sum  of  the  forcing  tides  (M2+N2)  by  68° 
(standard  deviation  error  of  five  degrees),  while  MS4  leads  the  sum  of  the 
forcing  tides  (M2+S2)  by  61°  (standard  deviation  error  of  four  degrees). 

These  phase  relationships  enhance  the  time  asymmetry  favoring  ebb.  The  con¬ 
sistency  of  the  phase  lead  in  all  overtides  and  compound  tides  throughout 
the  estuary,  and  for  different  periods  of  time,  suggests  it  is  a  property 
of  the  estuary  as  a  whole,  and  not  just  of  local  channel  characteristics. 

c)  Generation  of  the  forced  MSf,  or  fortnightly,  compound  tide:  The 
MSf  tide  consists  of  both  an  elementary  component  which  is  part  of  the  equi¬ 
librium  tide,  and  a  possible  forced  component  arising  from  the  interaction  of 
M2  and  S2  (S2-M2).  At  the  offshore  gage,  MSf  is  smaller  than  within  the 
estuary,  with  different  phase  relationships,  suggesting  that  the  more  ener¬ 
getic  MSf  found  in  the  estuary  itself  is  partly  a  forced  component. 

Although  other  longer-term,  less  energetic  constituents  with  similar  frequen¬ 
cies  exist,  this  study  concentrated  on  the  MSf  constituent  with  appropriate 
error  analysis,  as  it  is  the  largest  of  these  fortnightly  constituents.  Once 
inside  the  estuary,  MSf  grows  from  a  value  of  near  zero  to  between  five  and 
ten  centimeters.  There  is  no  consistent  pattern  of  growth  of  MSf  within  the 
estuary,  rather  its  amplitude  depends  in  part  on  atmospheric  conditions, 
suggesting  it  is  not  a  simple  compound  tide  (table  4).  Of  three  periods  of 
observation,  the  second  period  uniformly  had  the  largest  relative  amplitudes 
of  MSf  (compared  to  the  total  tide),  and  it  was  the  most  energetic  period 
from  a  storm  perspective.  The  first  and  third  periods  were  characterized 


by  lower  wind  stress.  As  an  illustration,  the  tidal  residual  which 
included  fortnightly  components  shows  relatively  low  amplitude  in  the  first 
period,  more  energetic  storm  super-elevation  during  the  second  period,  and 
relatively  less  energy  during  the  third  period  (fig.  10).  In  each  of 
these  residuals  a  clear  fortnightly  structure  is  apparent,  with  super¬ 
imposed  oscillations  of  various  frequencies.  Fourier  analysis  of  resi¬ 
duals  shows  the  fortnight] y  band  to  be  cuspate,  particularly  during 
periods  of  large  wind  stress  events,  illustrating  the  complexity  of  the 
interactions  among  tidal  and  ron-tidal  forcing. 

The  MSf  phase  provides  some  insight  into  the  nature  of  the  observed 
fortnightly  tide,  although  phase  information  is  more  complex  than  for  other 
compound  tides.  During  the  two  periods  of  lower  mean  wind  stress  the  forced 
fortnightly  tide  is  in  phase  with  the  forcing  tide  (S2-M2)  with  little 
scatter  about  the  mean  phase  difference  (table  4).  During  these  periods, 

MSf  is  fairly  small  (except  at  Nauset  Bay).  During  the  middle  period, 
which  included  a  large  storm  event,  MSf  phase  locking  is  still  excellent 
(values  of  140°,  141°,  and  141°  for  the  three  estimates),  although  MSr 
leads  the  forcing  by  a  significantly  different  amount  than  during  the 
other  two  periods.  This  anomalous  phase  relationship  is  related  to  the 
superelevation  corresponding  to  the  9  October  1982  storm  event,  which 
increased  water  levels  by  approximately  0.75m  within  the  estuary.  The 
inconsistency  in  phase  relationships,  and  apparent  dependence  on  meteoro¬ 
logical  events,  suggests  that  prediction  of  the  MSr  fortnightly  tide  Is 
tenuous,  with  consequent  poorly  forecasted  mean  water  level  deviations  of 
10  cm  or  so  on  a  fortnightly  cycle. 

The  phase  locking  between  MSr  and  the  forcing  components  results  in  a 
low  water  level  which  is  lower  during  neap  tides,  higher  during  spring 
tides  (fig.  11).  This  behavior  has  been  observed  previously  in  rivers 
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gure  10.  The  fortnightly  tide  (MS.)  is  both  an  astronomic  and  forced  con¬ 
stituent  (figure  10a,  top  panel)  and  can  contain  energy  from  storm 
surge  (figure  10b,  bottom  panel).  These  two  examples  show  the 
residual  tide  when  diurnal  and  higher  frequency  constituents  are 
removed . 
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(e.g.,  LeBlond,  1979),  and  has  been  ascribed  to  a  fortnightly  modulation 
of  the  frictional  forces  due  to  fortnightly  variations  in  the  velocity 
field.  A  similar  mechanism  may  be  responsible  for  the  observed  fortnightly 
relationship  in  Nauset  and  similar  estuaries,  although  the  response  of 
Nauset  is  complicated  by  its  response  to  wind  stress.  This  fortnightly 
component  is  under  further  investigation  by  the  authors. 

V.  IMPLICATIONS  FOR  SEDIMENT  TRANSPORT 

Tidal  asymmetries  generated  in  estuaries  through  non-linear  interac¬ 
tions,  and  represented  as  over  tides  and  compound  tides,  can  have  profound 
effects  on  sediment  transport  in  these  systems.  For  the  shelf  around  the 
British  Isles,  Pingree  and  Griffiths  (1979)  found  sand  transport  directions 
coinciding  with  asymmetries  in  the  tide  caused  by  the  M2,  M«  interaction. 
Boon  and  Byrne  (1981)  discussed  the  potential  importance  of  inlet/estuarine 
tidal  response  to  long-term  sedimentation  trends,  as  well  as  the  different 
types  of  asymmetry  present  in  shallow  inlet/estuarine  systems. 

For  Nauset  Inlet,  the  observed  asymmetry  favors  flood  transport  of  sedi¬ 
ment,  and  promotes  infilling  of  the  estuary.  This  trend  is  depicted  in  the 
phase  relationship  between  the  overtides  and  compound  tides  with  the  domi¬ 
nant  elementary  tidal  constituents,  where  the  higher  frequency  induced 
tides  always  lead  the  elementary  constituents  in  sea  surface.  This  phase 
lead  in  sea  surface  results  in  flood  flows  which  are  shorter  in  duration, 
and  consequently  stronger,  than  ebb  flows.  Since  sediment  transoort  is 
proportionaL  to  the  excess  shear  stress  raised  to  some  power,  the  higher 
flood  velocities  may  result  in  a  net  up-estuary  (flood)  transport  of 
material,  providing  the  drag  coefficient  remains  nearly  equal  over  the 
tidal  cycle.  In  the  distal  reaches  of  tidal  channels  at  Nauset,  the 
increase  in  magnitude  of  flood  flow  is  sufficient  that  the  transport  is 
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greater  despite  the  longer  duration  of  the  less  intense  ebb  flow.  The 
balance  between  increased  flow  velocity  and  lesser  duration  for  other 
estuaries  may  not  be  as  clear-cut,  so  each  asymmetrical  tide  must  be  eval¬ 
uated  given  its  peculiar  geometry  and  dynamical  characteristics.  Although 
this  study  has  emphasized  sea  surface  elevation  instead  of  velocity  as  an 
indicator  of  asymmetry,  the  flow  records  available  are  consistent  with  this 
flood  asymmetry,  and  verify  the  landward  transport  tendency. 

The  existence  of  a  landward  transport  of  coarse,  near-bed  material  is 
reflected  in  the  geometry  of  the  far  reaches  of  the  estuary.  In  the  north 
(Nauset  Bay,  station  8),  the  mean  water  depth  is  much  less  than  in  other 
parts  of  the  estuary.  Although  partly  a  result  of  barrier  overwash,  the 
geometry  of  the  flood-tide  deltas  shows  a  clear  up-estuary  transport.  Even 
more  clearly,  the  channels  leading  into  Town  Cove  in  the  southern  part  of 
the  estuary  are  shoaling,  with  a  well-defined  flood-tide  delta  encroaching 
on  the  deeper  waters  of  Town  Cove.  Local  measurements  of  the  flow  struc¬ 
ture  near  this  delta  demonstrate  the  dominance  of  up-estuary  transport. 
Unless  conditions  change  drastically,  the  estuary  will  continue  to  fill, 
with  the  lower  mean  water  depths  resulting  from  this  filling  process  accel¬ 
erating  the  rate  of  fill. 

On  a  time  scale  of  decades  or  a  century,  the  increase  in  sea-level 
along  this  region  may  also  affect  the  estuarine  response  to  tidal  forcing. 
Sea-level  rise  in  this  region  has  been  approximately  2.5  mm/yr  over  the 
past  40  years  (Aubrey  and  Emery,  1983).  Predictions  of  increased  rate  of 
sea-level  rise  over  the  next  century  due  to  the  global  increase  in  atmos¬ 
pheric  carbon  dioxide  will  modify  this  pattern  significantly,  perhaps  even 
in  as  short  a  time  as  a  few  decades.  If  sea-level  rise  exceeds  the  rate 
of  filling  of  the  estuary,  the  system  may  develop  a  more  linear  response 
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to  tidal  forcing.  The  decreasing  rate  of  infilling  could  lead  to  a  more 
stable  estuary  in  the  future.  These  trends  are  presently  under  study  to 
evaluate  the  effects  of  increased  sea-level  rise  on  estuarine  physics. 

Fine  suspended  sediment  transport  will  be  affected  by  the  nonlinear 
tidal  distortion  as  well.  In  particular,  for  suspended  particulate 
material  consisting  of  plant  matter  or  fine  inorganic  sediment,  an 
increased  tendency  exists  at  Nauset  for  up-estuary  transport,  a  result  of 
the  flood  asymmetry.  Larger  flood  currents  exert  higher  near-bed  shear 
stresses,  and  hence  entrain  larger  quantities  of  material.  Sufficiently 
small  or  light  material  can  be  transported  in  the  suspended  mode.  The 
tidal  asymmetry  has  strong  implications  for  other  aspects  of  mass  flux 
calculations  commonly  attempted  in  estuarine  systems.  Many  investigators 
(e.g..  Boon,  1978;  Valiela  et  al.,  1978;  Kjerfve  et  al.,  1981;  Ward,  1981) 
have  attempted  to  measure,  with  varying  degrees  of  success,  dissolved  and 
particulate  organic  and  inorganic  species  over  tidal  cycles  to  determine 
net  influx  or  efflux  of  this  material  from  the  estuary.  Because  of  the 
volume  and  complexity  of  sampling  required,  and  logistical  considerations, 
these  experiments  generally  last  30  hours  or  less.  To  perform  mass  bal¬ 
ances,  the  water  flux  must  be  balanced;  that  is,  the  net  volume  of  water 
flowing  from  the  estuary  must  be  equal  to  the  sum  of  freshwater  inflow  and 
evapo-transp i rat  ion  terms  (which  sum  can  be  near  zero).  Since  these  latter 
terms  are  difficult  to  quantify  independently,  a  tidal  cycle  of  12.4  or 
24.8  hours  is  often  imposed.  Results  commonly  show  a  net  influx  or  efflux 
of  water  over  the  tidal  period,  which  must  be  accounted  for  in  the  chemical 
balances.  Primary  contributors  to  this  net  water  flux  are  diurnal  inequal¬ 
ity,  meteorological  events  and  the  MSr  (fortnightly)  compound  tide.  MSf 
at  Nauset  can  have  a  10  cm  amplitude  (total  range  of  20  cm),  accounting 
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for  an  appreciable  difference  in  water  storage  within  the  estuary  system 
between  Spring  and  Neap  tides.  If  measurements  of  mass  flux  are  made 
during  periods  between  Spring  and  Neap  tides,  much  of  this  water  overburden 
will  be  transported  either  into  or  out  of  the  estuary,  resulting  in  poten¬ 
tially  large  net  water  fluxes.  MSf  can  also  affect  estimates  of  residual 
flows.  The  period  from  neap  to  spring  tides  will  correspond  to  a  period 
of  net  landward  flowing  residual  currents,  and  spring  to  neap  will  corres¬ 
pond  to  net  seaward  flowing  currents.  If  experiments  are  performed  over 
peak  Spring  or  peak  Neap  tides,  this  effect  will  be  reduced.  Storm  sur?c 
associated  with  wind  stress  events  can  also  dominate  the  net  flux  calcula¬ 
tions,  as  superelevations  of  1  m  or  more  are  not  uncommon  during  large 
storms.  The  timing  of  these  storms  can  severely  bias  mass  flux  calcula¬ 
tions.  In  summary,  lower  frequency  tidal  and  non-tidal  motions  limit  the 
accuracy  of  mass  balances  in  estuaries.  Since  these  lower  frequency  flows 
contribute  significant  water  fluxes  through  an  estuary,  they  must  be  pro¬ 
perly  sampled  to  obtain  an  accurate  estimate  of  mass  balances,  a  difficult 
criterion  to  satisfy. 

VI.  CONCLUSIONS 

Observations  of  sea  surface  elevation  and  horizontal  currents  in  a 
shallow  inlet/estuarine  system  demonstrate  the  nature  of  the  tidal  distor¬ 
tion  and  suggest  implications  for  sediment  transport.  The  major  findings 
supported  by  these  observations  include: 

i)  Overtides  and  compound  tides  derived  from  the  dominant  equilibrium 
tidal  constituents  offshore  can  be  important  in  estuaries  such  as  Nauset 
which  are  relatively  shallow,  with  large  tidal  amplitudes,  and  moderate 
tidal  flat  extent.  The  primary  interaction  is  with  the  M2  tide,  whose  over¬ 
tides  dominate  the  non-linear  signature  of  the  estuary.  Growth  of  the  M2 
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overtides  as  well  as  the  compound  tides  demonstrates  the  degree  of  non¬ 
linearity,  since  these  constituents  continue  to  grow  within  the  estuary  in 
spite  of  strong  frictional  dissipation.  Because  of  the  filtering  character¬ 
istics  of  the  estuary,  the  M4  tide  is  the  dominant  overtide,  since  dissi¬ 
pation  at  higher  frequencies  is  greater  than  at  lower  frequencies  (the 
estuary  is  a  low  pass  filter,  to  first  order). 

ii)  The  overtides  and  compound  tides  at  frequencies  higher  than  semi¬ 
diurnal  are  phase  locked  to  their  forcing  tides,  with  the  forced  tides 
consistently  leading  their  forcing  by  60-70°  in  phase.  This  phase  lead 
appears  to  be  a  characteristic  of  the  overall  geometry  and  dissipation 
characteristics  of  the  inlet/estuarine  system,  and  does  not  reflect  any 
offshore  phase  relationship.  The  relative  phase  is  maintained  throughout 
the  estuary,  even  as  the  amplitude  of  the  forced  tides  increases.  Reasons 
for  this  phase-locking  are  discussed  in  a  companion  paper  by  Speer  and 
Aubrey  (1984). 

iii)  The  sense  of  the  phase  locking  (with  forced  tides  leading  the 
forcing  in  sea  surface)  imparts  a  flood  dominance  to  Nauset  Inlet,  such 
that  flood  tides  are  shorter  but  more  intense  than  ebb  tides.  This  char¬ 
acteristic  can  be  inferred  from  the  sea  surface  records,  where  flood  tide 
is  shorter  and  ebb  tide  is  longer,  as  well  as  from  direct  observations  of 
currents  within  the  estuary.  Estuaries  with  different  phase  relationships 
can  have  ebb  dominance,  where  ebb  tides  are  shorter  but  more  intense  than 
flood  tides.  Watchapreague  Inlet,  Virginia,  an  example  of  an  ebb-dominated 
estuary  (sea  surface  M4  to  M2  phase  of  approximately  -140°  to  -160°;  Boon 
and  Byrne,  1981),  is  used  as  a  case  example  in  the  companion  theory  paper 
by  Speer  and  Aubrey  (1984). 
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iv)  A  fortnightly  compound  tide  is  found  consistently  within  the  estuary. 
The  MSf  constituent  is  partly  generated  within  the  estuary  by  interaction  of 
the  S2  and  M2  tides,  perhaps  modulated  by  frictional  effects,  as  suggested 

by  LeBlond  (1979)  for  riverine  flow.  This  10  cm  overtide  causes  lowest 
water  levels  at  Neap  tides  rather  than  at  Spring  tides,  since  MSr  is  in  phase 
with  its  forcing  (for  two  of  the  three  periods  of  observation).  The  ampli¬ 
tude  and  phase  relationships  are  consistent  throughout  the  estuary  for  any 
given  observation  period,  so  the  estuary  appears  to  respond  as  a  system. 
During  one  of  three  observation  periods,  a  large  storm  induced  a  0.75  meter 
superelevation  within  the  estuary,  with  a  decay  time  of  approximately  one 
week.  This  energetic  event  swamped  the  non-linear  component  of  MSr,  such 
that  the  storm  surge  dominated  estimates  of  MSf.  Even  though  estuarine 
response  was  uniform  at  the  two  working  tide  stations  within  the  estuary, 
phase  relationships  and  amplitudes  were  much  different  than  for  non-storm 
periods.  This  latter  observation  provided  estimates  of  MSf  which  were 
therefore  dynamically  unrelated  to  the  compound  (S2-M2)  tide. 

v)  From  a  sediment  transport  standpoint,  the  overtides  and  compound 
tides  at  frequencies  greater  than  semi-diurnal  cause  net  import  of  near¬ 
bed  material  at  Nauset.  More  intense  flood  flows  transport  more  sediment 
than  less  intense  ebb  flows  of  longer  duration  (aided  by  the  threshold 
criterion).  This  causes  filling  of  the  estuary,  as  evidenced  by  shoaling 
in  the  extreme  parts  of  the  estuary  (notably  Town  Cove).  Over  a  long-time 
frame,  decreased  depths  may  accentuate  the  filling  tendency,  although  the 
effects  of  increased  sea  level  rise  due  to  increased  global  atmospheric 
carbon  dioxide  may  reverse  this  trend.  For  suspended  particulate  trans¬ 
port,  the  asymmetry  experienced  by  the  Nauset  inlet/estuary  system  favors 


up-estuary  transport,  as  higher  flood  shear  stresses  (on  the  average)  will 
entrain  more  fine  sediment  and  light  particulate  material  for  subsequent 
transport . 

vi)  The  effects  of  the  MSf  fortnightly  tide  are  more  subtle,  as  MSf 
has  a  longer  period  and  does  not  directly  affect  the  tidal  asymmetry.  MSf 
makes  accurate  estimates  of  mass  balances  more  difficult,  as  the  superim¬ 
posed  low-frequency  flow  can  impart  a  significant  bias  to  water  flux  esti¬ 
mates  over  a  25  hour  period.  Mixing  is  substantially  influenced  as  well, 
since  MSf  can  account  for  a  large  volume  of  water  (given  its  20  cm  range). 
Estimates  of  residual  currents  must  consider  the  effects  of  these  lower  fre¬ 
quency  flows. 
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CHAPTER  IV 

A  STUDY  OF  NON-LINEAR  TIDAL  PROPAGATION  IN  SHALLOW 
INLET/ESTUARINE  SYSTEMS 
PART  II:  MODELING  TIDAL  ASYMMETRIES 
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Abs tract 

The  development  of  tidal  asymmetries  is  studied  via  numerical  integra¬ 
tion  of  the  one-dimensional  equations  for  channel  geometries  characteristic 
of  those  in  shallow  estuaries  located  on  the  U.S.  east  coast.  Channels 
without  tidal  flats  develop  a  time  asymmetry  characterized  by  a  longer 
falling  than  rising  tide.  This  behavior  is  enhanced  by  strong  friction 
and  large  channel  area  variability  over  a  tidal  cycle.  The  resulting  tidal 
currents  have  a  shorter,  intense  flood  and  a  longer,  weak  ebb  (flood  domi¬ 
nant).  Addition  of  tidal  flats  to  the  channels  can  produce  a  longer  rising 
tide  and  stronger  ebb  currents  (ebb  dominant),  if  the  area  of  tidal  flats 
is  large  enough  to  overcome  the  effects  of  time-variable  channel  geometry. 
Weaker  friction  with  flats  can  also  produce  this  asymmetry.  The  general 
features  of  tidal  channels  in  shallow  estuaries  can  be  represented  with 
simple  geometries.  Despite  the  physical  complexity  of  these  systems, 
essential  features  of  the  estuarine  tidal  response  can  be  recovered  from 
simple  models.  This  is  illustrated  by  modeling  idealized  representations 
of  tidal  channels  at  Nauset  Inlet,  MA  and  Wachapreague  Inlet,  V A. 
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1.  Introduction 

The  tidal  response  of  shallow  inlet/estuarine  systems  commonly  found 
on  the  U.S.  east  coast  is  frequently  characterized  by  the  development  of 
asymmetries  in  the  estuarine  tide.  The  asymmetries  take  the  form  of  un¬ 
equal  duration  and/or  unequal  magnitude  of  ebb  and  flood.  This  distortion 
of  the  estuarine  tide  can  be  represented  by  the  non-linear  growth  of  har¬ 
monics  of  the  astronomical  constituents.  A  particularly  important  type  of 
tidal  distortion  is  the  duration  inequality  in  the  rise  and  fall  of  the 
surface  tide.  Some  inlet/es tuarine  systems  (e.g.,  Murrels  Inlet,  SC,  Boon 
and  Byrne,  1981;  Nauset  Inlet,  MA,  Aubrey  and  Speer,  1983)  are  character¬ 
ized  by  longer  ebbs  than  floods  and  consequently  a  tendency  for  higher 
current  velocities  during  flood  (flood  dominant).  Other  systems  (e.g. 
Wachapreague  Inlet,  VA,  Boon  and  Byrne,  1981;  North  Inlet,  SC,  Nummedal 
and  Humphries,  1978;  Price  Inlet,  SC,  FitzGerald  and  Nummedal,  1983)  show 
the  opposite  duration  asymmetry  and  hence  a  tendency  for  stronger  ebb  cur¬ 
rents  (ebb  dominant).  The  different  types  of  asymmetry  have  important 
implications  for  estuarine  sediment  transport,  dispersal  of  water  column 
contaminants,  and  (over  long  time  scales)  es tuar ine/inlet  stability.  An 
important  question,  therefore,  is  what  processes  are  responsible  for  differ¬ 
ent  forms  of  tidal  asymmetry.  This  paper  describes  a  series  of  numerical 
experiments,  based  on  idealized  estuarine  geometries,  which  examine  the 
problem  of  estuarine  tidal  distortion. 

2.  Physical  Description  of  the  Problem 

Before  proceeding  to  the  numerical  solutions,  a  description  of  phys¬ 
ical  aspects  of  the  problem  is  presented.  For  representative  systems, 

(fig.  1.),  the  estuary  is  serviced  by  a  narrow  inlet  and  consists  of  a 
number  of  channels  generally  terminating  in  bays.  The  bays  may  consist 
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largely  of  tidal  flats  (fig.  la)  or  may  be  basins  deeper  than  the  connect¬ 
ing  tidal  channels  (fig.  lb).  Separating  the  tidal  channels  are  regions  of 
flats  exposed  at  low  tides  and  islands  of  marsh  grass  which  are  occasion¬ 
ally  overtopped  at  high  tide.  Typical  channel  depths,  h,  are  1-10  meters 
and  lengths,  L,  are  on  the  order  of  1-10  km.  Characteristic  values  for  the 
tidal  amplitude,  a,  suggest  that  finite  amplitude  effects  as  indicated  by 
the  ratio,  a/h,  can  be  important  in  these  systems.  Since  the  tidal  wave¬ 
length,  2ir/k,  is  much  greater  than  L,  characteristic  length  scales  are 
determined  by  estuarine  geometry.  The  channels  are  also  long  compared  to 
their  width,  b/L<<l,  and  the  horizontal  aspect  ratio,  h/b,  is  small.  The 
flow  may  therefore  be  termed  one-dimensional  although  channel  bends  visible 
in  Figure  1  will  produce  locally  two-dimensional  effects  (e.g.,  Aubrey  and 
Speer,  1984). 

Estuarine  systems  of  interest  to  this  study  are  also  characterized  by  a 
well-mixed  water  column  and  negligible  freshwater  inflow.  Mixing  is  pro¬ 
moted  by  the  strong  currents  (order  1  m/s)  often  found  in  estuarine  channel 
and  by  interaction  of  currents  with  a  rough  bed.  In  deeper  basins  at  chan¬ 
nel  terminations,  a  weak  stratification  develops  during  summer  months. 

The  offshore  astronomical  tide  spectrum  is  composed  of  a  number  of 
important  constituents  (fig.  2a).  A  complex  interaction  between  these 
constituents  is  possible  in  a  non-linear  system.  At  Nauset,  field  obser¬ 
vations  of  estuarine  tides  indicate  a  rich  spectrum  of  forced  motions 
(fig.  2b)  ranging  from  high  frequency  (i.e.,  greater  than  semi-diurnal)  to 
low  frequency  (e.g.,  fortnightly).  The  M2  constituent,  however,  as  the 
major  constituent  along  the  northeast  U.S.  coast  dominates  non-linear  pro¬ 
cesses  within  the  estuary.  This  domination  can  be  shown  from  theoretical 
arguments  (e.g.,  Gallagher  and  Munk,  1971)  as  well  as  field  observations 
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Figure  2a.  Periodogram  derived  from  tidal  harmonic  analysis  of  pressure 
sensor  located  in  10  m.  water  depth  outside  Nauset  Inlet. 
Subscripted  letters  refer  to  tidal  constituents  with  M2  the 
dominant  constituent.  Relative  size  of  the  ratios  M4/M2  and 
M6/M2  specifies  the  degree  of  non-linearity  of  the  tides.  In 
this  case  ratios  are  low  indicating  a  weakly  non-linear  tide. 
Ratios  are  for  amplitudes  not  energy. 
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Figure  2b.  Periodogram  of  tides  within  Nauset  estuary.  Note  the  growth  of 
high  frequency  (>2  cpd)  overtides  and  compound  tides.  Hit.  large 
M4/M2  ratio  indicates  strong  distortion  of  the  surface  tide. 
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(Aubrey  and  Speer,  1983).  The  harmonics  of  M2  are  larger  than  compound 
tides  produced  by  interaction  among  constituents. 

Field  observations  (e.g.,  Boone  and  Byrne,  1981;  Aubrey  and  Speer, 

1983)  and  theoretical  work  (e.g.,  Kreiss,  1957;  Dronkers,  1964;  Pingree 
and  Griffiths,  1979)  have  shown  that  the  interaction  of  M2  and  its  first 
harmonic,  M4 ,  explain  the  general  features  of  observed  tidal  asymmetries. 

The  type  of  tidal  distortion  (e.g.,  flood  dominant  versus  ebb  dominant) 
depends  on  the  relative  phasing  of  M4  to  M2.  Defining  the  M2  and  M4 
constituents  (sea  surface  or  velocity)  as 

Am  =  ai  cos(cot  -  0i  ) 

2 

Am  =  a2  cos(mt  -  02) 

4 

the  M2  -  M4  phase  can  be  written: 

<j>  =  201  -  0  z . 

The  magnitude  of  tidal  distortion  is  determined  by  the  rate  of  harmonic 
growth  and  is  indicated  by  the  ratio  of  M4  to  Mz  amplitudes: 

M4/M2  =  a2/ai 

An  important  problem  is  determining  the  estuarine  characteristics  which  con¬ 
tribute  to  non-linear  tidal  response.  This  study  investigates  the  mechanics 
of  M2  propagation  in  channels  with  dimensions  similar  to  those  described 
previously.  The  one-dimensional  equations  of  motion  (cross-sectionally 
averaged)  are  numerically  integrated  for  various  examples  of  idealized 
estuarine  channels.  Although  this  is  clearly  a  simplification  of  the 
actual  physical  problem,  the  above  discussion  shows  that  insight  can  be 
gained  from  this  approach. 

Boundary  forcing  for  the  problem  consists  of  the  M2  constituent.  Har¬ 
monics  present  in  the  solution  will,  therefore,  be  entirely  locally  gener¬ 
ated  within  and  confined  to  the  model  channel.  This  situation  is  compli- 
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cate  ’  in  the  case  of  actual  inlet/estuarine  systems  by  the  presence  of 

* 

shallow  features,  ebb  tide  deltas,  seaward  of  the  inlets.  Finite  amplitude 
effects  on  the  ebb  tide  delta  produce  harmonics  of  M2  which  enter  the  t 

estuary  through  the  inlet.  Field  evidence  (e.g.,  Boone  and  Byrne,  1981; 

Aubrey  and  Speer,  1983)  suggests  that  phase  and/or  amplitude  characteris¬ 
tics  of  M4  can  change  dramatically  through  the  inlet.  The  input  spectrum 
is  strongly  modified  by  the  estuarine  geometry.  Estuarine  response  solely 
to  M2  forcing  is,  therefore,  an  important  part  of  the  tidal  asymmetry 
problem.  In  principle,  the  modeling  can  include  an  input  spectrum 
consisting  of  the  other  important  astronomical  constituents  (e.g.,  Oi ,  Ki  , 

S2,  N2).  However,  the  most  general  features  of  the  estuarine  channel 
response  to  tidal  forcing  are  revealed  by  considering  M2.  Ter-diurnal  com¬ 
pound  constituents  (e.g.,  MK3,  MO3)  do  not  produce  any  asymmetry  and  quarter- 
diurnal  compound  constituents  (e.g.,  MSa ,  MN4 )  have  smaller  amplitudes  than 
M4. 

3.  Model  Formulation 

3.1.  Governing  Equations 

In  this  section,  we  develop  the  equations  of  motion  for  one-dimensional 
flow  in  a  well-mixed  water  column,  and  apply  them  to  cross-sections  of 
typical  estuarine  channels  (figs.  3  and  4).  A  trapezoidal  geometry  is  used 
to  represent  estuarine  channels  where  the  width  generally  contracts  with 
decreasing  depth.  We  consider  the  tidally  driven  problem  and  focus  on 
cases  where  sea  level  fluctuations,  a,  are  of  the  same  order  as  the  water 
depth,  h.  This  is  usually  the  case  in  these  estuaries  (e.g.,  Kjerfve, 


197b). 
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Examples  of  channel  cross-sections 


Figure  4.  Typical  estuarine  channel  cross-sections  from  Price  Inlet,  S.C. 
(taken  from  FitzGerald  and  Nummedal,  1983). 
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The  following  assumptions  are  made: 

1)  hydrostatic  pressure 

2)  level  sea  surface  across  channel 

3)  small  horizontal  aspect  ratio*  h/b<<l 
h)  long  narrow  estuary,  b/L<<l 

Assumption  one  simplifies  the  z-momentum  equation  to: 

p  =  pg(C-z)  (1) 

Assumption  two  simplifies  the  y-momentum  equation: 

C  =  C  (x,t)  (2) 

Clearly,  in  real  estuaries,  these  assumptions  will  be  violated  locally. 
Channels  with  large  bedforms  and  shallow  depths  may  depart  from  purely 
hydrostatic  conditions.  Strong  channel  curvature  (e.g.,  flow  around 
bends)  may  violate  assumption  two.  This  two-dimensional  effect  is  gener¬ 
ally  localized,  affecting  secondary  circulation,  shear  stress  and  sediment 
transport  patterns  in  the  vicinity  of  the  bend  (Smith  and  McLean,  1983). 
This  work  examines  the  response  to  tidal  forcing  of  the  entire  estuarine 
channel.  Assumpt ions  3  and  U  hold  for  many  estuaries  of  interest. 

The  model  dynamics  will  be  expressed  in  terms  of  cross-sectional  flux, 
U(x,t),  cross-sect ional ly  averaged  velocity,  U(x,t),  and  sea  surface  eleva¬ 
tion,  (^(x,t).  Defining  the  velocity  vector  q(u,v,w),  the  governing  equa¬ 
tions  are: 

9u  +  3v  +■  3w  =  0  ( 3 ) 

3x  3y  3z 

9n  +  9  u2  +  9  uv  +•  9  uw  =  - 1  9p  +  1 
9t  9x  9y  9z  p  9x  p 


3T£jt. 

9z 


9t? 

9y 


Equation  (3)  is  integrated  vertically  over  the  water  column  and  horizontally 
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b2  C 

across  channel  (i.e.,  X  Xdzdy)  to  yield 

-b  i  -h 

3A  +  3  U  =  0  (5) 

3t  3x 

A(x,t)  =  cross-sectional  area  of  channel 

b2  C 

U  =  X  Xudzdy  =  CJ»A 
-b  i  -h 

For  a  trapezoidal  channel, 

3A  =  b  8£ 

3t  9t 

b(x,z)  =  channel  width  at  depth,  z. 

Continuity  becomes 

3£  +  1  3_  U  =  0  (6) 

3t  b  3x 

The  same  integration  on  x-momentum  yields  after  considerable  algebra 
(Appendix  1): 

3U  +_3  U2  =  -gA  3£  -  Tb  p  (7) 

3t  3x  A  3x  p 

Tb  =  average  shear  stress  on  solid  boundaries 
P  =  wetted  perimeter 

The  derivation  of  equation  (7)  requires  assumptions  to  integrate  the 
advective  and  frictional  terms.  The  form  of  the  advective  term  in  (7) 
depends  on  small  cross-sectional  variability  in  velocity.  This,  in  turn, 
requires  gentle  changes  in  cross-sec t ion  along  channel.  A  momentum  factor 
(e.g.  Dronkers,  1964)  may  be  introduced  to  account  for  non-uniformity  in 
velocity  across  the  channel.  For  the  purposes  of  this  study,  a  momentum 
factor  was  not  used.  The  friction  term  in  (7)  includes  sidewall  friction 
by  considering  channel  sidewalls  in  the  wetted  perimeter.  The  problem  is 
closed  mathematically  by  formulating  friction  as  follows: 
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Tb  =  pf|u|u 

f  =  dimensionless  friction  factor 

In  the  formalism  of  cross-sectional  discharge,  U,  bottom  friction  takes  the 
form, 

xb  =  pf lulu  (8) 

A*A 

This  form  is  one  commonly  used  by  modelers  to  parameterize  the  inherently 
non-linear  effects  of  bottom  friction  (e.g.,  Owen,  1980;  Uncles,  1981). 

Other  investigators  use  the  Chezy/Manning  expression  (e.g.,  Pearson  and 
Winter,  1976,  1977;  Teubner  and  Noye,  1978),  although  these  are  quadratic 
forms  as  well.  A  number  of  shortcomings  exist  in  this  formulation  of  fric¬ 
tion.  The  quadratic  form  is  strictly  valid  for  steady  flows.  Additionally, 
if  surface  gravity  waves  are  present,  a  constant  friction  coefficient  does 
a  poor  job  of  representing  bottom  friction  (Grant  and  Madsen,  1979).  Sur¬ 
face  gravity  waves  are  potentially  important  near  the  inlet  and  may  consid¬ 
erably  enhance  the  tidal  bottom  friction  in  this  location.  However,  this 
is  a  complication  which  will  not  be  treated  in  the  model.  The  modeling 
will  examine  tidal  propagation  in  the  absence  of  waves. 

Controversy  exists  as  to  whether  a  quadratic  form  accurately  models 
the  generation  of  tidal  harmonics.  In  particular,  it  is  felt  that  odd  har¬ 
monics  of  the  fundamental  frequency  are  over-emphasized  (Prandle,  1980). 
Snyder  e_t  aj..  (1979)  found  that  a  combination  of  quadratic  and  linear  fric¬ 
tion  was  required  to  reproduce  M4  and  M6  observations  in  a  shallow  gulf. 
Other  investigators  have  had  success  modeling  the  M2,  M4  interaction  uti¬ 
lizing  a  quadratic  friction  (e.g.,  Pingree  and  Maddock,  1978;  Pingree  and 
Griffiths,  1979).  This  numerical  study  is  not  an  attempt  to  reproduce  a 
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set  of  field  observations  but  rather  to  examine  the  types  of  tidal  asym¬ 
metry  produced  by  different  estuarine  geometries.  A  quadratic  friction 
term  is  adequate  for  the  purposes  of  this  study. 

Equations  (6)  and  (7),  together  with  the  form  of  bottom  friction  speci 
fied  in  (8),  constitute  the  equations  solved  by  the  numerical  model.  Prin 
cipal  non-linear  effects  in  these  equations  enter  through  quadratic  fric¬ 
tion,  advection  of  momentum  and  through  tidal  interactions  with  estuarine 
geometry  in  the  continuity  equation.  By  varying  channel  geometry  and  fric 
tion  characteristics,  the  interplay  of  these  terms  in  distorting  the  tidal 
signal  is  examined. 

3.2.  Inclusion  of  Tidal  Flats 

An  important  physical  characteris tic  of  these  systems  is  the  presence 
of  tidal  flats  exposed  at  low  tide.  Tidal  flats  can  occupy  a  significant 
fraction  of  the  total  estuarine  surface  area  (e.g.,  FitzGerald  and  Numme- 
dal,  1983).  They  can  be  incised  with  shallow  channels  which  are  active 
during  only  a  portion  of  the  tidal  cycle.  This  type  of  geometry  is  parti¬ 
cularly  prevalent  around  the  inlet  mouth  where  they  are  termed  marginal 
flood  channels  (Hayes,  1975).  Obviously,  a  one-dimensional  model  cannot 
handle  the  presence  of  several  distinct  channels  across  a  channel  sec¬ 
tion.  Throughout  most  of  the  estuary,  however,  tidal  flats  serve  mainly 
to  store  water  as  opposed  to  transporting  momentum.  Bottom  friction  sub¬ 
stantially  weakens  the  velocity  as  water  moves  away  from  the  main  channel 
onto  the  flats.  This  suggests  approximating  the  effect  of  tidal  flats  by 
including  them  solely  in  the  continuity  equation.  This  idealization  will 
neglect  dissipation  of  the  tide  as  water  flows  across  the  flats,  as  well 
as  the  dissipative  effect  of  low  velocity  fluid  draining  into  the  ebbing 
channel  flow  from  tidal  flats  during  falling  tide.  Dronkers  (1964) 
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Figure  5.  Channel  geometry  with  tidal  flats  included. 
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suggests  a  procedure  for  including  lateral  inflow  (discussed  below).  Con¬ 
sideration  of  flats  solely  as  reservoirs  has  been  utilized  successfully 
before  in  model  applications  (e.g.,  Teubner  and  Noye,  1978;  Boon  and  Byrne, 
1981). 

A  model  estuarine  channel  is  composed  of  a  momentum  transporting  portion 
of  maximum  width,  b0,  and  a  storage  region  of  initial  width,  2bt  (fig.  5). 

The  slope  of  the  storage  region,  tan0,  is  equal  to  Ah/bt  where  Ah  repre¬ 
sents  the  difference  between  total  channel  depth  and  the  depth  at  width  b0. 

The  presence  of  the  tidal  flats  enters  through  the  continuity  term,  f 
9U/9x,  where  b  is  the  total  width  including  flats.  The  dissipative  effect 
produced  by  water  on  the  flats  entering  the  main  channel  is  modeled  as  follows. 

The  lateral  inflow  of  tidal  flat  water  into  the  main  channel  per  unit 
length  of  the  channel  is  equal  to  (br-bo)  9C/9t,  where  bT  is  the  total 
width  of  the  channel.  The  acceleration  of  this  water  in  the  main  channel 
unit  area  equals: 

U  (b-r-bo  )9£ 

9t 

or 

U  (b-r-bo )9£  (9) 

A  3t 

This  term  will  carry  a  positive  sign  on  the  right  hand  side  of  (7).  Using 
continuity,  we  can  rewrite  (9)  as 

(b-r-bo)  U  9U  (10) 

bT  A  3x 

showing  that  this  term  is  quadratic  in  transport.  It  has  the  effect  of 
acting  like  additional  friction.  Unless  bT>>b0  this  term  is  insignificant 
compared  to  the  advective  term.  For  most  cases,  the  modeling  ignores  this 
dissipative  effect.  The  preceeding  derivation  can  also  be  obtained  by  re¬ 
taining  lateral  inflow  terms  in  momentum  and  continuity  (Appendix  I). 
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3.3.  Boundary  Conditions  and  Initial  conditions 

Equations  (6)  and  (7)  constitute  an  initial-boundary  value  problem. 

Initial  conditions  are  taken  as: 

U(x,0)  =  0 
C(x,0)  =  0 

The  strong  frictional  effects  in  the  models  quickly  damp  out  initial  trans¬ 
ients,  making  these  convenient  initial  conditions.  Boundary  conditions 
include  a  condition  of  no  flux  into  solid  boundaries  and  a  specified  tidal 
forcing  at  the  channel  open  boundary, 

U(L, t )  =  0 
C(0,t)  =  a  cos  ut 

3.4.  Numerical  Scheme 

The  numerical  integration  of  equations  (6)  and  (7)  is  carried  out  with 
an  explicit  time-stepping  finite  difference  model  (e.g.,  Teubner  and  Noye, 
1978;  Prandle,  1980;  Uncles  and  Jordan,  1980).  Sea  surface  and  transport 
are  spatially  staggered  with  width  and  depth  defined  at  sea  surface  grid 
points  (fig.  3).  Time  derivatives  are  replaced  by  forward  time  steps  and 
spatial  derivatives  by  centered  differences.  The  formal  accuracy  of  the 
numerical  scheme  is  approximately  0(At,  Ax2).  An  explicit  finite-differ¬ 
ence  scheme  was  chosen  because  the  relative  simplicity  and  small  size  of 
the  model  domain  made  it  easy  to  implement  and  inexpensive  to  run.  Another 
popular  scheme  for  tidal  problems  is  the  harmonic  decomposition  method 
(e.g.,  Pearson  and  Winter,  1976,  1977;  Snyder  e_t  al.  ,  1979;  Le  Provost  et 
al . ,  1981;  Uncles,  1981).  This  method  recasts  the  time-dependent  problem 
(hyperbolic  equations)  into  a  time-independent  problem  (elliptic  equations). 
The  resulting  spatial  problem  is  solved  via  finite  differences  or  finite 
elements.  An  attractive  feature  of  this  method  is  that  one  solves  directly 
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for  the  amplitudes  and  phases  of  tidal  harmonics.  It  is  particularly  use¬ 
ful  for  large  2-D  problems  which  are  linear  or  weakly  non-linear  and  for 
modeling  a  number  of  different  constitutents .  Non-linear  effects  are 
included  either  by  analytical  expansions  for  non-linear  terms  (e.g.,  Le 
Provost  e_t  al.»  1976,  Snyder  et^  al . ,  1979)  or  by  iterative  approaches 
(e.g.,  Pearson  and  Winter,  1976).  Strongly  non-linear  problems  (such  as 
those  considered  here)  may  converge  slowly  or  not  at  all.  In  fact,  Pearson 
and  Winter  (1976,  1977)  found  it  necessary  to  include  an  artificial  fric¬ 
tion  term  to  insure  convergence  in  strongly  non-linear  coases.  For  this 
reason,  as  well  as  the  ease  of  including  non-linear  terms,  a  time-stepping 
model  was  employed. 

The  continuous  variables  A,  U  and  P  are  replaced  by  their  discrete 
analogues  over  the  space  and  time  domain  Xj=jAx  and  tn  =  nAt  for  j  =  0,N 
and  n  =  O,00.  The  remaining  variables,  C  and  b,  are  replaced  by  a  discrete 
set  over  Xj+,/2  and  tn  for  j  =  0,N  and  n  =  O,00.  Linear  interpolation 
defines  values  of  the  discrete  variables  for  all  x  and  t.  Defining  the 
following  operators: 


N" 

=  N(x j ,tn 

) 

AtN  " 

=  (wri  - 

Nr/At 

AxNj 

=  (Nr,  - 

Nr /Ax 

SxN" 

=  (Nj* I / 2 

-  Nr,/Z)/Ax 

SXN j , ! 

i  =  (wr,  - 

Nj- i )/ 2Ax, 

equations  (6)  and  (7)  are  replaced  by  their  discrete  analogues: 

At£j«-l/2  +  1  j  +  1  /  2  *  AxU  j  =  0  (11) 

b 

AtU1}  +  5,0?. ,  = 

-gA rxcr1  -  f iu"iu"«(p£ 


) 


(12) 
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and  D*J  =  (U?)2/A?. 

Mean  velocity,  U,  is  defined:  U"  =  Uj/A" 

The  boundary  and  initial  conditions  become: 

U"  =  0  j  =  N  (13) 

C"  =  a  cosu(nAt)  j  =  1/2 

(14) 


and 

C  0 

ZD 

=  u?  +  b?/2 

•Ax*At£ 1/2 

and 

Uj 

=  C5+i/2=  0 

j  =  0 ,N. 

The  stability  criterion  governing  the  linearized  numerical 

the  well-known  CFL  condition  (e.g.,  Roache,  1972): 

Ax  >  /ghmax 
At 


problem  is 


This  condition  states  that  information  cannot  be  propagated  faster  than  a 
"computational  information  speed"  defined  by  the  length  and  time  scale  of  the 
model.  The  form  of  the  advective  term  used  in  this  numerical  scheme  can  cause 
numerical  instability  under  certain  circumstances.  The  simple  advective 
equation 

8 U  +  c3u  -  0, 

3t  3x 


c  =  advection  speed 

centered  in  space  and  forward  in  time  (FTCS),  is  unstable.  The  instability 
in  this  equation  can  be  removed  by  one-sided  rather  than  space-centered 
differencing  (Roache,  1972).  It  is  also  possible  to  stabilize  the  advec¬ 
tive  term  in  the  FTCS  scheme  via  the  "angled  derivative"  approach  (e.g., 
Flather  and  Heaps,  1975;  Uncles  and  Jordan,  1980).  This  method  performs  a 
forward  and  backward  sweep  of  the  grid  utilizing  two  different  forms  of  the 
advective  term,  resulting  in  an  advective  term  centered  in  space  and  time. 
Although  the  present  study  does  not  employ  this  scheme,  the  dissipation  in 
the  momentum  equation  overcomes  any  tendency  towards  advective  instability. 
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The  procedure  in  model  runs  was  to  time-step  a  particular  experiment 
until  a  periodic  solution  was  obtained.  The  amplitude  and  phase  relation¬ 
ships  of  harmonics  were  determined  by  harmonic  analysis  of  the  solution. 

Sea  surface  and  mean  velocity  were  Fourier-analysed  according  to: 

C  =  <C>  +  £(A„cosnut  +  Bnsin  nut) 
n 

U  =  <U>  +  £(Cncos  nut  +  Dnsin  nut) 
n 

<  >  =  tidally-averaged  quantity 

n  =  1  :  M2 

etc . 

n  =  2  :  M4 

Additionally,  the  terms  in  the  equations  of  motion  were  harmonically  analysed 
to  examine  mechanisms  for  particular  cases  of  non-linear  tidal  distortion 
(e.g.,  Prandle,  1980). 

3.5.  Model  Comparison  with  Linear  Solution 

Although  the  general  form  of  the  numerical  scheme  described  above  has  been 
employed  successf ully,  it  was  tested  against  a  linear  solution  for  tidal 
propagation  in  a  channel.  The  approach  taken  was  to  solve  the  full  model 
equations  (with  a  modified  form  of  friction)  for  a  case  expected  to  be 
approximately  linear.  This  solution  was  then  compared  to  a  solution  of  the 
linearized  equations  which  are,  of  course,  only  valid  in  a  weakly  non¬ 
linear  case.  Such  a  case  occurs  when  a/h<<l  and  when  strongly  time-variable 
geometry  is  excluded.  The  channel  used  for  the  comparison  was  20  km  long 
with  constant  width  (200  m)  and  depth  (20  m)  forced  by  a  tide  with  ampli¬ 
tude,  0.5  m.  The  ratio,  a/h,  equals  0.025  and  cross-sectional  area  changes 
only  5%  over  a  tidal  cycle.  The  linearized  equations  governing  motion  in 
t-he  channel  are: 
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1  3U  +  3C  =  o 

b0  3x  3t 

3U  =  -gAD3£  -  fU 

3t  3x 

b0  =  channel  width  (constant) 

A0  =  mean  channel  area  (constant) 
f  =  linear  friction  factor  (t'1) 

The  solution  for  sea  surface  subject  to  boundary  conditions 

U(L, t )  =  0 
C(0,t)  =  H  coscot 

and  initial  conditions, 

U(x,0 )  =  C(x,0)  =  0, 

is  (Dronkers,  1964): 

C(x,t)  =  H  A  coscot  +  B  sincot  (18) 

cosh(so)  +  cos(sB) 

where, 

A  =  2^i(L)*,(x-L)  +  4>2(L)4>2(x-L)]] 

B  =  -2[jML)<Mx-L)  -  4>2(L)<j),(x-L)3 

<J>i(x)  =  cosh(lsor)cos(TsI3r) 

4>2(x)  =  sinh(lsCTr)sin(-2-sBr) 

s  =  2Lco  Co  =  (gh0  )  1  /z 

Co 

o  =Q-0.5  +  0.5(1  +  f  2)‘'2[]1/2 

13  =  £0.5  +  0.5(1  +  f 

co 

The  numerical  model  solves  equations  (6)  and  (7)  with  a  linear  friction 
for  various  grid  spacings  and  time  steps.  The  effect  of  different  grid 
spacings.  Ax,  on  the  numerical/analytical  comparison  was  examined  (Table 
l  and  fig.  6).  The  maximum  time  step  for  a  particular  Ax  was  determined 
by  the  CFL  condition.  In  addition,  the  size  of  the  time  step  was  varied 


(16) 

(17) 


for  one  particular  value  of  Ax.  Results  (Table  1)  show  that  differences 


TABLE  1 


Expt. 

AX1 

L 

At2 

T 

Al3 

(xlO"4 5) 

A24 

(xlO'4) 

FDT1 

.0125 

0.84 

3.24 

5.62 

FDT2 

.0250 

0.84 

4.18 

7.78 

FDT3 

.0500 

0.84 

6.72 

13.34 

FDT4 

.0500 

0.42 

6.54 

13.12 

FDT5 

.0500 

0.21 

6.48 

13.00 

Numerical  experiments 

.0357 

0.36-0.44 

1.  Grid  spacing  divided  by  channel  length. 

2.  Time  step  divided  by  time  scale  derived  from  the  Courant  condition 
(T  =  AX/Vgh). 

3.  Rms  error  over  tidal  cycle  in  sea  surface  for  numerical  solution 
compared  to  analytical  solution  normalized  by  tidal  amplitude  (50  cm) 
at  0.45  length  of  the  channel. 

4.  Same  as  3.  but  at  0.9  length  of  the  channel. 

5.  Values  of  space  and  time  steps  for  numerical  experiments  described 
later. 


SEA  SURFACE  (CM) 


ELAPSED  TIME 

ANALYTICAL/NUMERICAL  COMPARISON  (FDT2)  -  GRID  POINT  4 


ANALYTICAL/NUMERICAL  COMPARISON  (FDT2)  -  GRID  POINT  20 

Figure  6.  Comparison  of  linear  analytical  solution  and  numerical  solution 
for  example  FDT2  (see  table  1) .  The  comparison  is  for  sea  sur¬ 
face  at  two  grid  points  within  the  model  channel. 
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between  the  numerical  solution  and  analytical  solution  are  small  for  the 
range  of  Ax  and  At  values  examined  As  can  be  expected,  the  numerical/ 
analytical  differences  decrease  for  decreasing  Ax  (increased  spatial 
resolution  of  the  channel)  and,  to  a  lesser  extent,  for  decreasing  At. 

A  comparison  of  sea  surface  at  two  grid  points  for  one  of  the  examples 
(FDT2)  shows  good  visual  agreement  in  amplitude  and  phase  (fig.  6).  The 
model  runs  described  subsequently  have  spatial  and  temporal  resolution 
(Ax  =  250  m.  At  =  15  s.)  within  the  ranges  studied  in  this  section. 

4.  Numerical  Experiments 

4.1.  Introduction 

The  numerical  experiments  fall  into  two  general  classes: 

1)  Tidal  propagation  in  channels  of  constant  mean  width  and  depth, 
both  with  and  without  flats. 

2)  Tidal  propagation  in  channels  which  represent  features  of  the 
geometry  of  two  real  estuaries. 

Each  class  was  examined  for  the  sense  of  tidal  asymmetry,  and  for  the 
phasing  and  rate  of  growth  of  the  first  harmonic  (M«). 

The  first  class  of  numerical  experiments  is  similar  in  philosophy  to  a 
number  of  analytical  efforts  which  have  examined  aspects  of  tidal  propaga¬ 
tion  in  channels.  Kreiss  (1957)  noted  the  occurrence  of  tidal  i  urrent  time 
and  velocity  asymmetries  in  estuaries.  He  studied  the  case  of  flood  cur¬ 
rents  shorter  in  duration  but  reaching  greater  peak  velocity  than  ebb 
(flood  dominant),  and  explained  this  behavior  as  a  result  of  advectional 
non-linearities  and  linear  frictional  effects.  This  work  was  extended  by 
GaLIager  and  Munk  (1971)  who  demonstrated  the  effects  of  quadratic  friction 
and  an  incident  tide  with  many  constituents.  Quadratic  friction  introduced 
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odd  harmonics  into  the  solution  at  second  order  in  a/h.  Kabbaj  and  Le 
Provost  (1980)  utilized  a  complicated  analytic  decomposition  of  quadratic 
friction  to  include  more  accurately  the  effects  of  friction.  These  exam¬ 
ples  are  limited  to  cases  of  simple  constant  width  and  depth  rectangular 
channels,  and  the  ratio  a/h<<l.  In  fact,  as  previously  discussed,  natural 
estuarine  channels  tend  to  have  widths  which  vary  with  depth.  To  include 
this  feature,  as  well  as  cases  where  a/h  is  not  constrained  to  be  much  less 
than  one  (again,  a  common  feature  of  these  estuaries),  numerical  integra¬ 
tion  of  the  equations  is  necessary. 

The  second  class  of  experiments  adds  features  characteristic  of  natural 
estuarine  systems.  These  experiments  compare  two  inlet/estuarine  systems 
with  different  tidal  responses.  The  two  estuaries  are  Nauset  Inlet,  MA 
(fig.  lb)  and  Wachapreague  Inlet,  VA  (fig.  la).  Nauset  Inlet  is  character¬ 
ized  by  strong  harmonic  growth,  particularly  of  the  Mq  constituent,  and  an 
asymmetry  favoring  Longer  ebbs  than  flood  (e.g.,  Aubrey  and  Speer,  1983). 
Wachapreague  Inlet  has  much  weaker  harmonic  growth  and  an  asymmetry  favoring 
longer  floods  (e.g.,  Byrne  et  ai.,  1977;  Boon  and  Byrne,  1981).  Each  estuary 
is  modeled  by  an  idealized  channel  which  represents  its  basic  geometry. 

These  models  are  a  combination  of  the  features  described  in  the  first  group 
of  numerical  experiments.  The  trade-off  in  effects  of  various  geometries  on 
the  development  of  tidal  asymmetries  can  be  studied  in  these  cases. 

l*.2.  Parameters  Studied  in  Model  Runs 

Referring  to  Figure  5,  the  following  paramenters  define  the  geometry  of 
the  problem: 
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h  =  total  depth 

h2  =  mean  channel  depth 

hi  =  channel  depth  at  width  b0 

a  =  tidal  amplitude 

b0  =  width  of  momentum  transporting  channel 

Ab  =  difference  between  width  at  channel  top  and  bottom 

bt  =  width  of  tidal  flats  at  depth  h2 

f  =  friction  factor  (parameterizes  bottom  roughness) 

L  =  channel  length 

Non-dimensional  parameters  describing  the  geometry  of  model  runs  can  be 
formed  from  this  group.  In  most  cases  the  channel  length  is  held  constant 
at  7  km,  a  characteristic  length  of  these  estuaries.  The  importance  of 
non-linear  terms  in  the  equations  of  motion  is  indicated  by  the  ratio  a/h. 
Characteristic  values  of  this  parameter  range  from  0.1  to  0.5  for  the 
estuaries  of  interest.  An  "amplitude  dispersive"  behavior  in  shallow 
water  waves  with  a  cascade  of  energy  into  harmonics  of  the  fundamental 
frequency  is  anticipated  for  this  range  of  values  (e.g.,  Gallagher  and 
Munk,  1971). 

Variability  in  cross-sectional  area  and  channel  width  over  a  tidal 

cycle  is  described  by  the  ratios  a/h  and  Ab/b.  This  is  shown  by 

considering  expressions  for  cross-section  and  width  for  the  channel 

depicted  in  Figure  3.  The  cross-section.  A,  and  width,  b,  are  given  by: 

A  =  A0  +  bo*C(l  +1  Ab«£_)  (19) 

2  b0  h2 

A0  =  channel  area  at  mean  depth  h2 

b  =  bo(l  +  Ab_»£_)  ( 20 ) 

bo  h  2 

strong  modulation  of  both  channel  area  and  width  for  non-rectangular  geome¬ 
tries  is  possible  if  the  ratio  Ab/b0  is  large.  No  variability  occurs  in 
width,  of  course,  for  the  rectangular  case.  For  the  non-rectangular  case 
and  a  tidal  motion  given  by 
(,  =  acoswt, 
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cross-sectional  changes  will  contain  both  the  original  frequency  and  its 
first  harmonic  at  0(a/h).  In  fact,  for  a  tidal  motion  containing  many 
harmonics , 

C  =  £ancosnwt  a1>>a2>>an 

a  large  number  of  interactions  are  possible,  suggesting  that  time  variable 
geometry  may  be  an  important  source  of  harmonic  growth.  Modulation  in 
channel  geometry  will  have  an  important  effect  on  both  the  non-linear 
continuity  term  and  on  the  friction  term.  Because  cross-sectional  area 
and  channel  perimeter  change  at  different  rates  over  a  tidal  cycle,  fric¬ 
tion  effects  may  be  strongly  modulated  by  geometry.  The  effects  of  fric¬ 
tion  are  also  dependent  on  the  value  of  the  friction  factor  used  in  model 
runs.  Tidal  propagation  is  studied,  therefore,  for  a  range  of  friction 
factors . 

The  geometry  of  the  tidal  flats  is  described  by  the  slope  of  the 
storage  banks, 

tanG  =  Ah/bt  where  Ah  =  h2-hi 

As  mentioned  earlier,  the  influence  of  tidal  flats  is  introduced  through 
continuity.  Also  the  parameter  (Af/Ac)  comparing  the  fractional  change  in 
the  storage  area  versus  momentum  transporting  portion  of  the  channel  describes 
the  channel  geometry  with  tidal  flats: 


Other  quantities  of  interest  such  as  potential  tidal  prism  and  mean  tidal 
current  can  be  obtained  from  the  geometric  parameters.  The  potential  tidal 
prism  is  the  volume  of  water  between  high  and  low  water  the  estuary  would 
admit  in  the  absence  of  tidal  decay  (or  amplification).  The  reduction  of 
potential  tidal  prism  could  be  examined  instead  of  the  decay  in  tidal 
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range.  In  terms  of  tidal  dynamics,  these  two  points  of  view  are  equiva¬ 
lent.  The  mean  current  is  defined  as  the  potential  prism  divided  by  the 
product  of  channel  cross-sectional  area  times  one-half  the  tidal  period 
(Swift  and  Brown,  1983).  This  quantity  is  defined  assuming  the  absence  of 
prism  reductions  (tidal  decay)  and  tidal  asymmetry. 

In  the  following  sections,  the  parameters  are  varied  systematically  to 
examine  their  effects  on  tidal  distortion  in  the  model  channel.  The  role 
of  channel  geometry  in  producing  distorted  tides  is  explored  by  varying  the 
ratios  db/b0  and  a/h.  Models  with  the  same  geometry  are  run  with  differ¬ 
ent  values  of  the  friction  factor.  Finally  the  effects  of  variable  tidal 
flat  area  relative  to  channel  area  are  studied. 

A. 3.  Constant  Mean  Width/Depth  Problems 
4.3.1.  Physical  Considerations 

Two  types  of  geometry  are  considered  in  this  section:  channels  with 
and  without  tidal  flats.  Harmonic  growth  depends  on  the  balance  between 
energy  extracted  from  the  fundamental  frequency  by  non-linear  friction, 
advection  and  continuity  and  the  energy  lost  by  friction.  These  problems 
have  both  a  linearly  decreasing  potential  prism  and  mean  velocity  down- 
channel  because  of  the  no-flow  interior  boundary  condition  (fig.  7).  For 
these  geometries,  friction  and  advection  of  momentum  will  be  strongest  near 
the  forcing  boundary.  Non-linear  advection  which  scales  with  the  channel 
length,  L,  will  be  relatively  unimportant  as  an  M<*  generation  mechanism 
compared  to  friction  in  the  model  runs.  In  contrast  to  friction  and  advec¬ 
tion,  non-linear  continuity  will  show  only  a  slow  decline  (if  any)  down- 
channel.  Near  the  interior  boundary  little  change  occurs  in  the  growth  of 
harmonics  in  sea  surface  and  harmonics  in  velocity  decrease  to  zero.  Model 
runs  reach  an  "equilibrium"  phase  relationship  between  M2  and  M4  in  both 
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sea  surface  and  velocity  (or  transport)  which  is  determined  by  the  geometry 
and  frictional  characteristics  of  the  channel.  A  representative  model  run 
(fig.  8)  shows  rapid  adjustment  in  sea  surface  M2-M4  phase  4><r»  near  the 
forcing  boundary  followed  by  more  gradual  changes  near  the  interior  boundary. 
The  ratio  of  sea  surface  M4  to  M2  amplitude,  (M4/M2)<,  also  increases 
rapidly  near  the  forcing  boundary  before  leveling  off  near  the  interior 
boundary. 

Differences  exist  between  the  model  behavior  described  above  and  that 
of  natural  estuarine  channels.  Cross-sectional  area  can  vary  widely  along 
real  channels  producing  large  velocities  well  within  the  estuary.  As  a 
result,  friction  and  advection  will  not  necessarily  show  the  monotonic 
downchannel  decline  exhibited  by  the  models.  The  advective  term,  although 
relatively  unimportant  in  the  models,  can  become  locally  significant  in 
natural  channels  due  to  rapid  changes  in  geometry  (e.g.,  flow  around  bends, 
sudden  channel  expansions).  It  can,  in  fact,  dominate  over  friction  as  an 
M4  generation  mechanism  in  deeper  channels  with  different  geometries  than 
those  under  consideration.  Pingree  and  Maddock  (1978)  have  shown,  for  ex¬ 
ample,  that  the  advective  term  is  more  important  that  friction  for  M4 
generation  in  the  English  channel. 

4.3.2.  Numerical  Experiments:  channels  without  tidal  flats 

These  experiments  are  organized  in  terms  of  increasingly  time  variable 
channel  cross-sectional  area  as  shown  by  the  parameter,  a/h  •  Ab/b0  (Table  2). 
This  non-dimensional  parameter  is  strictly  geometrical,  indicating  how 
variable  the  channel  cross-section  area  is  over  a  tidal  cycle.  Since 
channel  area  enters  into  dynamical  terms  and  the  non-linear  continuity 
term,  a  large  value  of  this  parameter  can  indicate  the  possibility  of  a 


M4/M2  M2-W4  PHASE  (DEO) 
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TABLE  2 


Expt . 

a/h 

Ab/bo 

a/h»Ab/b0 

f 

1 

M4/M2 

2 

m4/m2 

3 

4 

4>u 

%AH 

10B 

.22 

0 

.0200 

.005 

.016 

92 

75 

0 

10 

.32 

0 

— 

.0200 

.019 

.065 

81 

50 

0 

10A 

.32 

0 

— 

.0020 

.002 

.010 

117 

84 

0 

10C 

.39 

0 

— 

.0200 

.019 

.065 

79 

55 

0 

10D 

.46 

0 

— 

.0200 

.027 

.092 

74 

48 

0 

22 

.32 

.20 

.06 

.0200 

.021 

.073 

76 

42 

1 

24 

.32 

.27 

.09 

.0200 

.022 

.076 

75 

38 

1 

24A 

.32 

.27 

.09 

.0100 

.013 

.046 

85 

59 

0 

24B 

.32 

.48 

.15 

.0200 

.026 

.093 

72 

26 

4 

13G 

.20 

.75 

.15 

.0200 

.009 

.031 

85 

56 

0 

13B 

.26 

.75 

.20 

.0200 

.021 

.074 

78 

31 

3 

13 

.32 

.75 

.24 

.0200 

.037 

.124 

68 

11 

10 

13E 

.32 

.75 

.24 

.0100 

.026 

.092 

75 

29 

3 

13F 

.32 

.75 

.24 

.0020 

.008 

.026 

101 

73 

0 

13A 

.32 

.87 

.28 

.0200 

.044 

.140 

65 

4 

16 

13C 

.39 

.75 

.29 

.0200 

.040 

.132 

70 

17 

8 

13D 

.32 

.94 

.30 

.0200 

.050 

.150 

65 

2 

19 

131 

.32 

.94 

.30 

.0020 

.014 

.047 

93 

62 

0 

13J 

.43 

.75 

.32 

.0200 

.056 

.168 

64 

3 

7 

1.  Sea  surface  amplitude  ratio  at  grid  point  4. 

2.  Sea  surface  amplitude  ratio  at  grid  point  26. 

3.  Sea  surface  M2-M4  relative  phase  at  grid  point  26. 

4.  Velocity  M2-M4  relative  phase  at  grid  point  26. 

5.  %  decline  in  tidal  range  at  grid  point  28. 
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strongly  non-linear  tidal  response.  It  is  not,  however,  a  fundamental 
measure  of  the  non-linearity  of  the  long  wave  (or  shallow  water)  equations. 
This  is  provided  by  the  ratio  of  wave  particle  velocity,  u,  to  wave  phase 
velocity,  c.  For  the  case  of  long  waves,  this  ratio  equals  a/h,  the  non- 
dimensional  wave  amplitude.  The  condition  of  linearity  for  long  waves 
requires  that  u/c<<l.  If  the  ratio  is  small,  a  large  value  of  Ab/b0  will 
have  no  impact  on  the  non-linear  development  of  the  tide.  Because  a/h  is 
large  for  these  examples,  channel  geometry  will  influence  the  tidal 
response. 

All  these  model  runs  produced  time  asymmetries  characterized  by  a  longer 
falling  than  rising  tide.  The  maximum  time  asymmetry  of  this  type  for  a 
given  ratio,  (M4/M2)<-,  is  produced  by  an  M2-M4  relative  phase,  4>o  of  +  90°. 
As  can  be  seen  in  Table  2,  <j>^  ranges  approximately  30°  on  either  side  of  90°. 
The  sense  of  time  asymmetry  is  determined  by  the  M2-M4  phase  and  the  magni¬ 
tude,  by  the  M4/M2  ratio.  Results  from  the  experiments  show  (M4/M2)<r 
depends  strongly  on  geometry  and  friction  characteristics  of  the  channel. 

In  general,  more  time-variable  geometry  and  stronger  friction  produce  larger 
M4/M2  ratios.  The  increase  in  this  ratio  is  dominated  by  growth  in  M4  as 
opposed  to  decrease  in  M2. 

These  numerical  examples  also  produce  little  decline  in  tidal  range 
except  for  the  cases  of  strongly  time-variable  channel  cross-section  area 
(e.g.,  experiments  13,  13A,  13D).  This  behavior  is  expected  considering 
the  linear  problem  (fig.  9).  The  linear  solution  shows  friction  in  a 
rectangular  channel  is  not  effective  in  reducing  tidal  amplitude  over 
short  channel  length  scales.  Modeling  the  large  amplitude  decay  often 
seen  in  natural  systems,  therefore,  requires  a  reasonable  representation 


of  channel  geometry. 
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Mean  velocity  exhibits  more  variability  in  M2-M4  phase  than  does  sea 
surface.  If  M4  is  in  phase  with  M2  in  velocity,  ($u  -  0°),  the  case  of 
short,  strong  floods  versus  long,  slow  ebbs  occurs  (e.g.,  experiments  13A, 
13D).  From  a  strict  continuity  argument,  this  situation  results  when 
falling  tide  exceeds  rising  tide  in  duration.  The  experiments  approxi¬ 
mating  this  particular  phase  relationship  are  characterized  by  strong 
friction  and  time  variable  channel  geometry.  In  general,  the  experiments 
summarized  in  Table  2  show  a  tendency  towards  higher  velocity  floods  than 
ebbs  moving  downchannel . 

The  tidal  harmonic  signature  in  velocity  is  more  complicated  than  sea 
surface  and  hence  is  not  summarized  in  Table  2.  All  cases  showed  greater 
harmonic  growth  in  velocity  than  in  sea  surface  (i.e.,  larger  M4/M2  etc. 
ratios).  Experiment  13  (fig.  10)  shows  that  whereas  M4  steadily  grows  in 
amplitude  for  sea  surface,  it  first  increases  and  then  declines  in  mean 
velocity.  Initially,  non-linear  processes  pump  more  energy  into  this  har¬ 
monic  than  is  removed  by  friction.  Eventually,  frictional  effects  and  the 
constraint  on  velocity  imposed  by  the  interior  boundary  condition  reduce 
its  amplitude.  The  increase  in  (M4/Mz)u  is  initially  a  result  of  both  decay 
of  M2  and  growth  of  M4 .  Note  also  that  higher  harmonics  in  velocity  (e.g.. 
Mg,  Mg)  decline  before  M4  begins  losing  energy. 

The  effect  of  channel  geometry  on  sea  surface  response  is  clarified  by 
three  experiments  (fig.  11)  shown  in  order  of  increasingly  trapezoidal  geo¬ 
metry  (i.e.,  larger  cross-sectional  variability  over  a  tidal  cycle).  Each 
experiment  has  the  same  friction  factor  and  tidal  forcing  (a/h  =  0.32). 

Sea  surface  becomes  increasingly  non-linear  as  channel  cross-sectional  area 
becomes  more  time-variable.  Both  frictional  effects  and  flow-geometry 
interactions  produce  this  reult.  Cross-sectional  area  contains  a  second 


LINEAR  SOLUTION 
FOR  A  RECTANGULAR  CHANNEL 


S2  =  0.4  41-  2  4  ^ 

Figure  9.  Linear  analytic  solution  for  tidal  propagation  in  a  rectangular 
channel  (equation  18).  The  vertical  axis  represents  sea  surface 
amplitude  at  interior  boundary  normalized  by  the  amplitude  at 
the  open  boundary.  S2  is  a  non-dimensional  channel  length  (=2kL) 
Sj  is  a  non-dimensional  friction  measure  (=fT/2).  The  channel 
length  scale  for  the  model  runs  is  marked  as  a  dotted  line. 


log  (cm/s)2  log  cm 
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EXPERIMENT  13 
SEA  SURFACE 


GRID  POINT 


VELOCITY 


Figure  10.  Periodogram  derived  from  harmonic  analysis  of  numerical  solution 
for  model  run  13.  Top  panel  is  sea  surface  at.  three  grid  points, 

bottom  is  velocity  at  the  same  points.  The  numbers  along  the 
horizontal  axis  refer  to  M2  ar  d  its  harmonics  (i.e.  1=M2;  2=^; 
3=M6;  4=Mb). 
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harmonic  term  (equation  19)  and  width  a  first  harmonic  when  the  channel 
geometry  is  a  trapezoid.  The  magnitude  of  these  terms  increases  as  Ab/b0 
increases.  Also,  velocity  tends  to  be  higher  in  the  trapezoidal  cases  as 
a  result  of  the  way  tidal  prism  and  channel  cross-sectional  area  combine. 

A  calculation  of  the  potential  prism  and  mean  velocity  (as  defined  previ¬ 
ously)  for  two  of  the  cases  (Table  3)  illustrates  the  point.  Higher  veloci¬ 
ties  lead  to  larger  values  of  friction  (assuming  f  =  constant)  and  hence 
more  harmonic  growth.  The  phase  relationship  between  M2  and  M4  in  each 
case  leads  to  falling  tide  exceeding  rising  tide  in  duration  with  the  asym¬ 
metry  increasing  downchannel.  The  distortion  of  the  tide  at  the  model's 
interior  boundary  for  three  examples  (fig.  12)  shows  the  following:  (i) 
the  duration  of  rising  tide  decreases  with  increasing  Ab/b0;  [For  experi¬ 
ment  10,  rising  tide  lasts  5  hr  40  m  (fall:  6  hr  45  m),  whereas  for  experi¬ 
ment  13D,  the  duration  of  rise  is  5  hr  10  m  (fall:  7  h  15  m).  The  latter 
value  is  comparable  to  time  asymmetries  observed  in  the  field  (e.j.,  Nauset 
Inlet,  Aubrey  and  Speer,  1983a)];  (ii)  the  reduction  in  tidal  range  asso¬ 
ciated  with  increasingly  trapezoidal  geometry;  and  (iii)  the  tendency 
towards  a  small  set  up  ( tidal ly-averaged )  in  sea  surface  with  increasingly 
trapezoidal  geometry.  This  non-linear  behavior  has  been  noted  before  in 
channel  models  by  Gallagher  and  Munk  (1971)  and  is  in  accordance  with  the 
concept  of  a  radiation  stress  (Longuet-Higgins  and  Stewart,  1964).  The 
greater  the  dissipation  of  the  tide,  the  larger  the  set-up  in  mean  sea 
level  at  the  channel  interior  boundary. 

Since  channel  cross-sectional  variability  over  a  tidal  cycle  is  also  a 
function  of  a/h,  increasing  this  ratio  will  lead  to  a  more  distorted  tide. 
Sea  surface  line  spectra  for  three  different  values  of  a/h  (fig.  13)  in 
cases  where  Ab/b  was  held  constant  show  enhanced  non-linear  transfers  to 
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TABLE  3  • 

POTENTIAL  PRISM/MEAN  VELOCITY 


Expt . 


potential  prism1 
(x  10 6m3) 


grid  point 


mean  velocity2 
(m/s ) 


M2  velocity3 

(m/s) 


Experiments  without  flats 


10 

1.26 

4 

0.17 

0.27 

12 

0.12 

0.18 

13D 

2.02 

4 

0.32 

0.40 

12 

0.22 

0.26 

Experiments 

with  flats 

11 

4.06 

4 

0.63 

0.58 

18 

0.23 

0.23 

14 

4.72 

4 

0.63 

0.55 

18 

0.26 

0.22 

1.  Potential  prism  (defined  in  text). 

2.  Mean  velocity  (defined  in  text). 

3.  M2  velocity  calculated  from  model  run 
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Figure  11.  Periodogram  of  sea  surface  for  three  model 
runs  at  two  different  grid  points.  These 

experiments  illustrate  the  effect  of  channel 
geometry  on  sea  surface  tidal  distortion. 


\\ 
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Figure  12.  Distortion  of  the  sea  surface  tide  near  model's  interior 
boundary  for  the  three  experiments  shown  in  figure  11. 
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all  harmonics  of  M2  as  a/h  increases.  Experiment  13J  had  a  larger  tidal 
amplitude  than  the  other  two  cases  in  which  the  ratio  was  increased  by 
decreasing  the  mean  depth.  Similar  to  the  geometry  parameter  Ab/b0,  larger 
values  of  a/h  lead  to  stronger  harmonics  in  cross-sectional  area  and  width 
(equations  19,  20).  As  mentioned  previously,  it  also  indicates  the  impor¬ 
tance  of  non-linear  terms  in  the  equations  of  motion.  Finally,  changing 
values  of  the  parameter  result  in  changes  in  mean  velocity.  This  leads  to 
increases  in  the  magnitude  of  the  advective  term  and  particularly  in  fric¬ 
tion.  The  first  two  examples  admit  approximately  the  same  potential  prism; 
however,  mean  cross-sectional  areas  decrease  from  13G  to  13.  Mean  channel 
velocities  are  larger  in  13  than  13G,  and  more  tidal  distortion  occurs  in 
the  former  than  the  latter  case.  Experiment  13J  admits  a  larger  prism 
than  the  others  and  also  has  larger  velocities.  Comparing  the  results  of 
experiments  24B  and  13G  (Table  2)  emphasizes  that  a/h  is  the  fundamental 
indicator  of  non-linearity  in  estuarine  tidal  response.  The  two  experiments 
have  approximately  the  same  value  of  the  geometry  parameter,  a/h  •  Ab/b0,  how¬ 
ever  experiment  24B  has  the  larger  a/h  ratio  and  the  more  non-linear  tide. 
Channel  geometry  effects  are  only  important  if  a/h  is  sufficiently  large. 

Mosc  of  the  experiments  were  run  with  friction  factors  of  0(10"2). 
Numerical  models  of  shelf  and  shallow  sea  circulation  commonly  use  a  coeffi¬ 
cient  of  0 ( 1 0 _  3 )  (e.g.  Flather  and  Heaps,  1975;  Beardsley  and  Haidvogel, 

1981),  although  this  value  is  too  low  if  surface  gravity  waves  are  present 
(Grant  and  Madsen,  1979).  Estimates  of  friction  factors  in  shallow  estuaries 
with  strong  tidal  flows,  large  bedforms  and  high  sediment  transport  rates 
are  of  the  order  10~2  (e.g.,  Swift  and  Brown,  1983;  Aubrey  and  Speer,  in 
prep.).  Results  of  experiments  with  3  values  of  f  ( . 002 , . 02 , . 04 ) ,  and  the 
same  geometry  and  forcing  (fig.  14)  show  the  effect  of  variable  f.  Friction 
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factors  of  1-2x10’ 3  (commonly  stated  values)  are  too  low  for  these  systems 
when  strong  tidal  flows  are  present.  This  level  of  friction  does  not  pro¬ 
duce  any  tidal  decay  and  also  causes  negligible  harmonic  growth. 

It  has  been  stated  that  quadratic  friction  produces  only  odd  harmonics 
(i.e.,  Ms)  and  that  even  harmonics  (i.e.,  M4 )  are  produced  by  non-linear 
continuity  and  advection  (e.g.,  Hamilton,  1977).  Analytical  models  employ¬ 
ing  perturbation  schemes  to  include  quadratic  friction  find  odd  harmonics 
generated  by  friction  to  second  order  in  non-dimensional  wave  amplitude 
when  a  time-invariant  friction  factor  is  used  (e.g.,  Gallagher  and  Munk, 
1971;  Kabbaj  and  Le  Provost,  1980): 

U  =  ecoscot  (21a) 

U|U|  =  T-rr  e2  {coswt  +  |cos3cot  - . }  (21b) 

where  e  represents  a  non-dimensional  measure  of  current  amplitude. 

Higher  order  interactions  (i.e.,  0(e3))  between  the  fundamental  frequency 
and  its  harmonics  are  not  calculated  in  analytical  models.  In  fact,  as 
Gallagher  and  Munk  (1971)  point  out,  the  even  harmonics  can  be  modified 
considerably  by  higher  order  frictional  interactions.  The  numerical 
results  presented  here  and  those  of  other  numerical  experiments  (e.g., 
Pingree  and  Maddock,  1978;  Prandle,  1980)  show  the  importance  of  friction 
in  generating  both  even  and  odd  harmonics.  This  is  important  because  even 
harmonics  (M4)  produce  the  time  asymmetries.  Distortion  of  the  tide  at  the 
model's  interior  boundary  as  a  function  of  friction  (fig.  15)  shows  a  pat¬ 
tern  similar  to  that  produced  by  increasing  the  magnitude  of  the  geometry 
parameter.  Although  increasing  friction  from  .02  to  .04  produces  a  more 
non-linear  tidal  response,  it  also  Leads  to  a  decrease  in  the  energy  of  the 
harmonics  and  a  stronger  attenuation  of  Mz  (fig.  14). 


Figure  14.  Periodogram  of  sea  surface  for  three  model 
runs  illustrating  the  effect  of  increasing 

the  friction  factor  on  tidal  distortion. 
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As  stated  previously,  the  mean  velocity  signal  is  more  non-linear  than 

sea  surface  in  all  these  experiments  (e.g.,  fig.  10).  Friction  and  advec- 

tion  act  directly  on  velocity  to  produce  harmonics  of  M2.  Conceptually, 

the  harmonic  growth  is  transferred  to  sea  surface  via  the  term,  1.  3A«U3. 

b  3x 

Since  b  has  a  constant  term  and  oscillates  dominantly  at  the  M2  frequency 
(equation  20),  the  effect  of  this  term  is  to  decrease  the  contribution  of 
the  higher  frequencies  in  sea  surface  relative  to  mean  velocity  or  transport. 
This  is  illustrated  as  follows: 

1  3  A«U  =  1  3  U 
b  3x  b  3x 

with, 

b  =  b0  +  bi  cos(ut  -  4>i  ) 

U  =  Ui  cos(wt  -  <j>2)  +  U2  cos(2iot  -  4>3) 

1  3U  =  Uix+Uzx  +  .  .  . 

b  3X  bo+bi 

=  l.(Ui  X+U2  x  +  .  .  .  .  )  ( 1  — b  i  + .  •  .  .) 

bo  bo 

Collecting  terms  of  frequency  u, 

Jjx  (22a) 

b0 

and  terms  of  frequency  2u, 

Uz_x  -  Ui,  •  b!  (22b) 

bo  b! 

values  from  experiment  13  show  that  whereas  the  ratio  of  U2/Ui  equals  0.27, 
the  ratio  of  the  2nd  and  1st  harmonics  of  the  continuity  term  (equations 
22b,  22a)  is  approximately  0.13.  The  pattern  in  experiment  13  is  similar  to 
the  others  and  shows  that  tidal  distortion  in  velocity  is  stronger  than  sea 


surface . 
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Similar  to  sea  surface,  increased  channel  cross-section  variability, 
non-dimensional  tidal  forcing  and  friction  lead  to  an  increasingly  distort¬ 
ed  mean  velocity  signal.  As  the  friction  factor  increases,  $u  shifts  towards 
zero  degrees  (table  2).  This  phase  relationship  also  requires  strongly  time- 
variable  channel  cross-sectional  area  (table  2).  Increasing  the  geometry 
parameter,  a/h  •  Ab/b0,  while  keeping  friction  constant  increases  the  magni¬ 
tude  asymmetry  and  the  duration  of  ebb  (fig.  16).  Note  that  velocity  is 
non-dimensionalized  by  the  maximum  velocity  developed  in  the  model  run.  The 
experiments  with  a  phase,  <f>u,  approximately  zero  are  also  the  ones  charac¬ 
terized  by  the  strongest  harmonic  growth  in  sea  surface  and  longest  falling 
tides . 

Three  experiments  (24,13,13D)  were  studied  in  more  detail  to  examine 
the  development  of  tidal  asymmetries  (tables  4,  5).  These  experiments 
focus  on  the  effects  of  geometry  (with  strong  friction)  in  producing  tidal 
asymmetry.  The  tidal  momentum  balance  is  dominated  by  pressure  gradient 
and  bottom  friction  (fig.  17).  Advection  of  momentum  is  not  an  important 
component  of  the  balance  and  the  temporal  derivative  is  only  important  dur¬ 
ing  the  turn  of  the  tide.  This  particular  balance  is  characteritic  of  all 
the  numerical  experiments. 

Non-linear  terms  in  the  equations  of  motion  were  Fourier  analysed  to 
provide  phase  and  amplitude  information.  The  form  of  the  terms  analysed 
in  this  section  are  (Appendix  I): 

f|ll|U  •  (s)  (friction) 

9  U2  (advection) 

9x 

1  3_  A»U  (continuity) 

b  9x 


2  3 


O  N  4  • 

*111 


UI0013A  (OSTTVIWON 


velocity  developed  in  that  model 
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TABLE  4 
FRICTION  TERM 


Expt.  Grid  km 

Location 


M2 (m/s2)  M4/M2  Relative  Phase 

M2-M4(°) 


24 

2 

0.5 

8.35 

X 

10'4 

.294 

90 

4 

1.0 

7.07 

X 

10'4 

.294 

83 

12 

3.0 

3.16 

X 

10'4 

.311 

66 

18 

4.5 

1.25 

X 

nr4 

.321 

61 

13 

2 

0.5 

1.76 

X 

10‘3 

.220 

104 

4 

1.0 

1.43 

X 

10~3 

.206 

84 

12 

3.0 

6.29 

X 

10'4 

.264 

41 

18 

4.5 

2.52 

X 

10'4 

.299 

33 

13D 

2 

0.5 

2.4 

X 

10'3 

.213 

117 

4 

1.0 

1.84 

X 

10'3 

.178 

85 

12 

3.0 

7.85 

X 

10’4 

.278 

28 

18 

4.5 

3.21 

X 

10-4 

.333 

20 
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TABLE  5 

CONTINUITY  TERM 


Expt . 

Grid 

Location 

km 

M2 (m/ s 2 ) 

m„/m2 

Relative 

M2-M4 

24 

2 

0.5 

1.26  x  10"4 

.056 

102 

4 

1.0 

1.25  x  10'4 

.075 

121 

12 

3.0 

1.24  x  10~4 

.139 

144 

18 

4.5 

1.24  x  10‘4 

.162 

149 

13 

2 

0.5 

1.24  x  10'4 

.138 

97 

4 

1.0 

1.20  x  10'4 

.165 

110 

12 

3.0 

1.14  x  10'4 

.256 

130 

18 

4.5 

1.13  x  10'4 

.293 

135 

13D 

2 

0.5 

1.22  x  10'4 

.172 

99 

4 

1.0 

1.15  x-  10*4 

.204 

113 

12 

3.0 

1.05  x  10‘4 

.304 

130 

18 

4.5 

1.04  x  10'4 

.354 

134 
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Harmonic  analysis  of  the  advective  term  confirms  that  it  is  small  compared 
to  friction  and  not  an  important  driving  mechanism  for  M<*  in  these  models. 

In  general,  the  M4  acceleration  from  friction  is  at  least  an  order  of 
magnitude  larger  than  that  from  advection.  This  conclusion  is  not  affected 
by  changing  geometry.  Although  the  magnitude  of  the  advective  acceleration 
increases  from  24  to  13D  in  response  to  the  larger  velocities  in  13D,  the 
relative  importance  of  advection  to  friction  does  not  increase.  Tidal  dis¬ 
tortion  is  controlled,  therefore,  by  friction  and  by  tidal  f low/geometry 
interactions . 

The  friction  term  is  strongest  in  experiment  13D  and  weakest  in  24 
(table  4).  This  is  a  result  of  both  the  higher  velocities  in  13D  and  the 
smaller  hydraulic  radius,  R  (=A/P).  For  the  three  experiments  at  grid 
point  4,  friction  is  larger  during  flood  than  ebb,  reflecting  the  velocity 
asymmetry  developed  in  these  models  (fig.  18).  The  hydraulic  radius 
(fig.  18)  varies  more  for  experiment  24  than  for  13  and  13D.  Large 
variability  in  R  arises  in  cases  where  channel  width  change  with  depth  is 
small  (i.e.,  channels  are  more  "rectangular”).  Despite  lower  velocities, 
experiment  24  (less  variable  geometry)  has  a  larger  frictional  M4/M2  ratio 
than  13  and  13D.  This  behavior  largely  results  from  the  greater  variability 
in  R  in  experiment  24. 

The  first  order  frictional  M4  generation  represents  interactions  between 
tidal  fluctuations  in  channel  geometry  at  the  Mz  frequency  and  frictional 
retardation  of  the  Mz  tidal  flow.  Between  the  time  of  high  water  slack  and 
peak  ebb,  friction  is  weaker  than  for  the  case  of  a  constant  mean  R.  This 
is  especially  true  when  R  shows  large  variability,  and  is  equivalent  to  a 
downchannel  acceleration.  Similarly,  between  low  water  slack  and  peak 
flood,  friction  will  be  stronger  than  for  the  case  of  constant  mean  R. 


IMS  (R) 


Figure  18.  Tidal  frictional  oscillation  for  three  experiments  at  grr 
point  4  (1  km).  The  form  of  friction  represented  in  this 
plot  is  (f/u/u) /R.  Variation  in  the  hydraulic  radius,  R, 

is  shown  on  the  plot. 
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This  is  again  equivalent  to  a  down-channel  acceleration.  Channel  geometry 
(entering  through  R)  modulates  friction  to  produce  an  M4  oscillation  with 
peak  down-channel  accelerations  between  low  water  and  maximum  flood,  and 
between  high  water  and  peak  ebb.  This  qualitative  explanation  describes 
only  the  first  order  friction-geometry  interaction.  Other  harmonics  in 
velocity  and  geometry  interact  in  a  complex  manner  not  readily  described. 
These  other  interactions  produce  the  different  M2-M4  phase  relationships 
observed  in  friction  in  the  experiments.  All  three  have  similar  phase 
relationships  (~+90°  -  100°)  near  the  forcing  boundary  (table  4). 

Further  downchannel,  experiments  13  and  13D  tend  towards  +20°,  whereas 
experiment  24  reaches  +60°  (table  4). 

Fourier  analysis  of  the  continuity  term  reveals  a  different  behavior 
than  the  friction  term.  The  magnitude  of  this  term  shows  little  decline 
downchannel,  and  in  all  three  cases,  harmonics  of  Mz  monotonically  increase 
downchannel  (table  5).  In  contrast  to  the  friction  term,  experiments  13D 
and  13  are  more  non-linear  than  24.  This  reflects  the  greater  variability 
in  channel  cross-sectional  area  over  a  tidal  cycle  and  greater  tidal  dissi¬ 
pation  in  the  former  two  experiments. 

The  phase  relationships  in  sea  surface  and  velocity  (4><,  4»u)  observed 
in  the  model  results  are  determined  primarily  by  friction  and  continuity. 

In  the  examples  described  above  and  the  experiments  listed  in  Table  2,  both 
terms  drive  the  sea  surface  phase  relationship  to  approximately  +90°  or 
less.  Strong  friction  and  time-variable  geometry  produce  a  sea  surface 
phase  of  ~65°  and  velocity  phase  of  ~0°  as  well  as  strong  harmonic  growth. 
Reducing  friction  (e.g.  experiments  13  vs.  13E,  13F)  for  the  same  variable 
geometry  changes  phase  relationships  dramatically.  Similarly,  making  the 
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geometry  less  variable  while  keeping  the  same  friction  (e.g.  experiments 
13  to  10)  causes  large  phase  changes.  Both  changes  also  have  the  effect 
of  reducing  harmonic  growth. 

4.3.3.  Channels  with  Tidal  Flats 

These  experiments  include  storage  regions  (tidal  flats)  with  some  of 
the  channel  geometries  described  in  4.3.2  (table  6).  Compared  to  the 
previous  experiments,  an  important  effect  of  the  storage  regions  is  to 
increase  the  potential  prism.  Purely  from  continuity  considerations, 
therefore,  larger  mean  velocities  are  anticipated  (table  3).  Frictional 
effects  will  be  enhanced  resulting  in  greater  dissipation  of  the  tide.  As 
in  the  previous  examples,  the  potential  prism  and  mean  velocity  linearly 
decrease  downchannel.  Again,  friction  and  advection  are  strongest  near 
the  forcing  boundary  whereas  the  non-linear  continuity  term  shows  little 
downchannel  decline. 

The  way  in  which  tidal  flats  are  introduced  causes  a  discontinuity  in 
the  channel  geometry  with  a  sharp  break  between  transported  and  stored 
fluid.  This  is  reflected  in  velocity  records  where  a  distinct  "break"  in 
the  record  occurs  as  tidal  flats  first  become  covered  (fig.  19).  The 
hypsometric  curve  of  two  experiments  (fig.  20)  also  emphasizes  the  break 
in  the  topography.  Real  estuaries,  of  course,  have  smoother  hypsometric 
curves.  The  ratio,  Ar/Ac,  (table  6),  indicates  the  change  in  tidal  flats 
cross-sectional  area  versus  channel  area  for  a  change  in  sea  level.  Three 
geometries  from  the  first  group  of  experiments  (10,  13,  13A)  were  utilized 
in  these  experiments  to  assess  the  effects  of  tidal  flats.  The  experi¬ 
mental  results  both  in  terms  of  phase  and  harmonic  growth  show  a  different 


TABLE  6 
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Velocity  M2-M4  relative  phase  at  grid  point  26. 
%  decline  in  tidal  range  at  grid  point  28. 


Figure  19.  Velocity  at  grid  point  4  from  model  run  11.  The  arrow 
refers  to  the  point  at  which  tidal  flats  become  flooded 

at  this  location.  Velocity  is  normalized  by  maximum 
velocity  developed  in  the  model. 


Figure  20.  Hypsometric  curves  for  model  runs  13  and  17.  The  sharp 

break  in  the  curve  for  17  indicates  the  presence  of  tidal 
flats. 


-167- 


pattern  than  the  previous  examples.  In  general,  patterns  of  harmonic 
growth  are  more  complex  for  these  cases  (the  complexities  will  be  discussed 
later)  but  certain  general  trends  can  be  derived. 

In  some  of  the  experiments,  the  sea  surface  phase,  <j>«r,  is  greater  than 
+180°.  This  phase  relationship  produces  a  longer  rising  tide  than  falling 
tide,  and  contrasts  with  the  previous  experiments.  The  channel  character¬ 
istics  include  lower  time-variability  in  channel  cross-section  and  large 
values  of  Ar/Ac.  Additionally,  the  results  of  experiments  14B  and  14C 
versus  14  suggest  that  weaker  values  of  friction  can  produce  a  longer  rise  even 
in  cases  of  more  variable  geometry.  In  general,  trapezoidal  channels  with 
relatively  small  tidal  flats  show  phases,  4>o  between  +100°  and  +180°  indi¬ 
cating  a  longer  falling  tide.  All  these  experiments  show  a  greater  adjust¬ 
ment  in  <j>c  than  the  previous  experiments.  In  every  example,  falling  tide 
exceeds  rising  tide  in  duration  near  the  forcing  boundary  (<J>C  <180°).  Where 
tidal  flats  are  important,  this  inequality  is  reversed  moving  downchannel 
(compare  fig.  8  with  fig.  21).  Thus,  large  storage  regions  adjacent  to  the 
channel  can  produce  a  duration  asymmetry  characterized  by  a  longer  rising 
tide. 

Harmonic  growth  and  particularly  the  relative  phase  of  M4  in  sea  surface 
depends  on  the  size  of  tidal  flats,  time  variability  of  channel  geometry  and 
the  level  of  friction.  The  first  six  experiments  listed  in  Table  6  involve 
rectangular  (i.e.  "weakly"  time-variable)  channels  with  different  tidal 
flat  areas  and  levels  of  friction.  The  geometry  is  similar  to  experiment 
10  (table  2)  which  has  a  longer  falling  than  rising  tide  (<^  -  +80°).  The 
addition  of  tidal  flats  to  this  channel  geometry  results  in  a  longer  rising 
tide  (i.e.  $«r  between  180°  and  360°).  As  the  channel  geometry  becomes  more 
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Figure  21.  Adjustment  in  sea  surface  phase  for  model  run  1 

(compare  with  figure. 8).  Below  the  dashed  line,  time 
asymmetry  is  characterized  by  a  longer  falling  tide; 

above  the  Line  by  a  longer  rising  tide. 
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time-variable  (through  larger  a/h),  the  phase  approaches  180°.  For  these 
experiments,  a  longer  falling  tide  at  he  model's  interior  boundary  does 
not  occur.  Adjustment  in  the  phase  of  sea  surface,  4>o  is  much  larger  in 
these  experiments  than  the  ones  in  Table  2  (fig.  21).  Because  of  this,  the 
sense  of  time  asymmetry  changes  moving  downchannel. 

A  comparison  of  the  surface  tide  at  grid  point  26  (fig.  22)  show  the 
different  type  of  distortion  present  in  experiment  10  (no  tidal  flats)  and 
experiments  11,  11A  (tidal  flats).  For  11  and  11A,  the  rise  of  water 
level  slows  as  soon  as  the  flats  begin  to  flood  (at  -10  cm).  This  causes 
longer  lags  in  high  tide  downchannel  compared  to  channels  without  flats. 

As  water  floods  the  tidal  flats,  pressure  gradients  driving  the  flow  relax 
and  the  rate  of  increase  in  velocity  slows  or  can  even  temporarily  reverse. 
This  imparts  a  distinctive  "kink"  in  model  velocity  records  (figs.  19,23). 
This  feature  has  been  observed  in  velocity  and  pressure  gradient  records  at 
Nauset  (Aubrey  and  Speer,  1983),  where  it  corresponds  to  tidal  flat  flood¬ 
ing.  The  velocity  signal  for  grid  ponts  downchannel  (e.g.,  grid  point  18) 
includes  not  only  the  signature  of  tidal  flat  flooding  but  also  reversals 
in  the  rate  of  velocity  increase  just  after  high  and  low  water.  These 
reversals  are  caused  by  pressure  gradient  fluctuations  driving  the  weak 
flows  in  the  model's  interior.  Note  that  these  features  are  not  present 
in  high  velocity  regions  rf  the  grid  (fig.  19).  They  do  not  appear  to  be 
directly  related  to  tidal  flat  flooding  and  uncovering,  however,  the  fluc¬ 
tuations  become  more  pronounced  with  increasing  tidal  flat  width.  If  the 
excursions  in  tidal  flat  width  become  too  large  for  a  given  channel 
b  =  bQ  +  bt  coscot 

and 

bt  >>  b0 , 

numerical  instabilities  develop  at  these  points  in  the  velocity  records 


Figure  23.  Velocity  at  grid  points  12  and  18  for  model  run  11. 
The  velocities  are  normalized  as  before. 
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and  eventually  contaminate  the  whole  solution.  A  limit  in  the  size  of 
flats  which  can  be  included  with  a  particular  channel  exists  beyond  which 
a  non-physical  solution  is  obtained. 

The  spectral  structure  of  sea  surface  at  grid  point  26  (6.5  km)  and 
velocity  at  grid  point  12  (3  km)  for  experiments  10,  11,  11A  shows  that 
the  tidal  flats  experiments  have  a  more  non-linear  response  (fig.  2A). 

The  effect  of  the  large  M4  component  in  the  tidal  flats  experiments  is 
apparent  in  the  time  asymmetries  of  both  sea  surface  and  velocity.  A 
strong  third  harmonic  term  is  forced  in  both  sea  surface  and  velocity  as 
well.  The  phasing  of  the  third  harmonic  in  velocity  causes  the  flattening 
of  the  curve  (fig.  23)  during  early  flood.  This  corresponds  in  part  to  the 
onset  of  flooding  of  the  tidal  flats.  The  M6  harmonic  is  more  nearly  in 
phase  with  sea  surface  M2  and  its  effect  is  not  as  noticeable  on  sea  surface 
records. 

The  strong  harmonic  growth  is  anticipated  by  the  larger  velocities 
associated  with  experiments  11  and  11A.  However,  generation  cf  third  har¬ 
monics  is  more  vigorous  for  tidal  flats  experiments  than  for  cases  of 
channels  without  flats  and  with  roughly  equal  velocities  (e.g.,  experiment 
13J,  table  2).  Examination  of  the  harmonic  structure  of  the  primary  non¬ 
linear  forcing  terms,  friction  and  continuity  for  experiments  10  and  11A 
reveals  the  difference  in  non-linear  response  (table  7).  With  strong 
dissipation  of  both  mean  velocity  and  sea  surface  near  the  forcing  bound¬ 
ary,  the  continuity  term  in  11A  drives  both  a  strong  M4  and  Ms  constituent. 
This  term  is  a  weak  source  of  non-linear  forcing  in  experiment  10  (compare 
also  with  table  5).  Conversely,  the  friction  term  in  experiment  10,  while 
having  a  smaller  absolute  magnitude,  nevertheless  is  a  more  non-linear  term 
as  measured  by  the  M4/M2  ratio. 


log  (cm/s) 


Figure  24.  Periodogram  of  sea  surface  and  velocity  for 

three  models  showing  the  effect  of  increasing 
tidal  flats  area  on  the  spectrum  of  the  estuarine 
tidal  response. 
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TABLE  7 

CONTINUITY  TERM 


Expt . 

Grid 

Location 

km 

M2(m/s2) 

m4/m2 

Relative 
M2-M4 ( 

Phase 

°) 

11A 

2 

0.5 

2.46  x  10"4 

.  446 

.120 

90 

10 

2.5 

2.32  x  10‘4 

.330 

.114 

87 

15 

3.75 

2.31  x  10'4 

.302 

.122 

80 

10 

2 

0.5 

1.26  x  10~4 

.014 

.005 

139 

10 

2.5 

1.26  x  10'4 

.057 

.019 

159 

15 

3.75 

1.26  x  10'4 

.100 

.033 

165 

Expt . 

Grid 

Location 

km 

FRICTION  TERM 

M2 (m/sz ) 

m4/m2 

Relative 

m2-m4( 

Phase 

°) 

11A 

2 

1.8  x  10'3 

.244 

.099 

200 

9.8  x  10'4 

.209 

.097 

208 

15 

6.1  x  10'4 

.200 

.107 

211 

10 

2 

6.0  x  10'4 

.346 

.072 

87 

4.3  x  10'4 

.350 

.071 

79 

15 

7.3  x  10'4 

.359 

.077 

71 
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The  next  group  of  experiments  listed  in  Table  6  introduces  the  effects 
of  more  time-variable  channel  geometry  with  variable  friction  and  tidal 
flat  area.  The  channel  geometry  is  essentially  the  same  as  that  used  in 
experiment  13  (table  2).  Tidal  asymmetries  developed  in  these  experiments 
reflect  the  conflicting  effects  of  channel  geometry  in  continuity  and 
momentum  tending  to  drive  a  longer  falling  tide  and  the  presence  of  tidal 
flats  in  continuity  driving  a  longer  rising  tide.  As  a  result,  a  variety 
of  channel  responses  including  both  types  of  time  asymmetry  are  found. 

Again,  an  important  point  about  these  model  runs  is  that  sea  surface  phase 
adjustment  between  M2  and  M4  can  be  substantial  (fig.  25).  This  is  particu¬ 
larly  true  of  the  experiments  which  involve  a  mixture  of  large  tidal  flats 
and  strongly  time-variable  geometry  (e.g.,  experiment  14).  The  important 
point  is  that  the  phase  relationship  does  not  level  off  quickly  as  in  the 
experiments  without  tidal  flats  (fig.  8).  The  phase  relationship  of  both 
sea  surface  and  mean  velocity  shown  in  Table  6  is  taken  from  near  the 
interior  model  boundary  and  indicates  the  type  of  tidal  asymmetry  the  model 
is  developing. 

In  general,  most  of  these  experiments  have  a  longer  falling  tide  (4>c 
<180°),  although  the  actual  phase  is  quite  different  from  experiment  13. 

The  size  of  tidal  flat  area  has  an  important  effect  on  the  type  of  tidal 
asymmetry.  Sea  surface  at  grid  point  26  for  experiment  14  shows  little 
time  asymmetry,  although  considerable  distortion  is  present  (fig.  26).  If 
the  amount  of  tidal  flat  area  is  increased,  the  example  of  a  longer  rising 
tide  is  found  (experiment  14A).  Conversely,  as  the  area  of  tidal  flats 
declines,  the  behavior  of  the  model  becomes  more  similar  to  experiment  13 
as  one  would  expect  (compare  experiments  17B  and  13).  A  comparison  of  sea 
surface  distortion  for  experiments  14A,  14  and  17B  illustrates  graphically 
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the  effect  of  combining  tidal  flats  with  time  variable  channel  geometry 
(fig.  26).  The  spectral  structure  of  sea  surface  for  these  examples 
reveals  the  strong  second  and  third  harmonic  growth  driven  by  cases  with 
extensive  tidal  flats  (fig.  27). 

The  rate  of  harmonic  growth  in  sea  surface  changes  as  the  area  of  tidal 
flats  decreases  (fig.  28).  The  case  of  extensive  flats  (experiment  14A) 
shows  strong  growth  of  both  M4  and  Ms.  As  the  area  decreases  (e.g.,  experi¬ 
ment  17)  harmonic  growth  is  sharply  reduced  and,  in  fact,  tends  to  decline 
downchannel  (fig.  28).  The  exact  reasons  for  this  behavior  are  not  clear, 
except  that  the  competing  non-linear  effects  of  tidal  flat  flooding  and 
time-varying  channel  geometry  must  be  nearly  in  balance.  With  continuing 
decrease  in  tidal  flat  area  (experiment  17B),  the  behavior  begins  to 
resemble  experiment  13  both  in  terms  of  harmonic  growth  and  phase  (fig.  28 
and  table  6).  The  pattern  of  M4  growth  and  decline  exhibited  in  experiment 
17  changes  to  one  of  steady  M4  growth.  Non-linear  effects  associated  with 
time-variable  channel  geometry  dominate  those  associated  with  the  presence  of 
tidal  flats  for  this  experiment. 

The  final  three  experiments  listed  in  Table  6  also  emphasize  the  inter¬ 
action  of  channel  geometry  and  tidal  flats.  The  larger  ratio  of  Ab/b0  indi¬ 
cates  greater  time-variability  of  channel  geometry.  The  channel  geometry  is 
the  same  as  experiment  13A  (table  2).  For  these  examples,  the  effects  of 
variable  channel  geometry  dominate  the  tidal  response  even  if  the  tidal  flat 
area  is  large  (experiment  23A).  The  result  is  that  all  three  examples  have 
a  tidal  asymmetry  characterized  by  a  longer  falling  tide.  Experiments  23A 
and  23  also  display  the  previously  observed  pattern  of  initial  M4  growth  and 
then  decay  representing  conflicting  effects  from  tidal  flats  and  variable 


channel  geometry. 


SEA  SURFACE 
GRID  POINT  26  (6.5km) 
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Figure  27.  Periodogram  of  the  surface  tide  for  the  grid  point  and 
three  experiments  shown  in  figure  26. 


CD 

^  Is-  Is- 


-181- 


The  magnitude  of  friction  within  a  channel  affects  its  tidal  asymmetry, 
as  seen  by  comparing  experiments  14,  14B,  14C  where  the  friction  factor  was 
reduced  while  the  geometry  remained  the  same.  Harmonic  growth  in  sea  sur¬ 
face  increases  slightly  with  increasing  friction  and  the  amount  of  tidal 
dissipation  increases  dramatically.  The  M2-M4  phase  of  sea  surface  shifts 
from  a  nearly  symmetrical  tide  in  experiment  14  to  a  longer  rising  tide  as 
the  level  of  friction  decreases  (table  6). 

4.4.  A  Comparison  of  Two  Different  Estuarine  Geometries 

The  problems  studied  in  the  previous  two  sections  identified  important 
estuarine  channel  characteristics  involved  in  different  types  of  asymmetry 
and  different  rates  of  harmonic  growth.  The  models  considered  in  this 
section  involve  two  different  geometries  based  on  estuaries  found  along 
the  U.S.  east  coast.  The  first  model  idealizes  the  case  of  the  main  tidal 
channel  at  Nauset  (MA).  The  second  is  approximately  representative  of  an 
estuarine  channel  at  Wachapreague  inlet/estuary  (VA).  These  numerical 
experiments  do  not  attempt  to  model  predictively  or  hindcast  conditions  at 
the  estuaries  in  question.  The  complexity  of  these  features  cannot  be 
reproduced  adequately  by  one-dimensional  models.  Instead,  the  models  are 
employed  as  diagnostic  tools  in  the  same  spirit  as  the  previous  sections. 
Actual  channel  depths  and  widths  are  idealized  and  not  reproduced  exactly. 
Model  results  are  discussed  in  terms  of  the  form  of  tidal  asymmetry  and 
not  in  terms  of  exact  matches  with  field  data.  The  problems  studied  here 
combine  elements  of  the  models  studied  in  the  previous  two  sections. 

Nauset  and  Wachapreague  represent  two  different  types  of  estuarine  geo¬ 
metry  with  different  tidal  responses.  Nauset  Inlet  is  characterized  by 
shallow  (<3  m)  tidal  channels  which  lead  to  an  unstable  inlet  with  flood 
and  ebb  tide  delta  deposits  (fig.  lb).  Tidal  flats  along  the  main  southern 
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drainage  are  located  near  the  inlet  mouth  and  consist  of  flood  tide  delta 
and  barrier  overwash  deposits.  The  total  extent  of  tidal  flats  is  moderate 
in  this  system.  Nauset  is  modeled  by  a  channel  of  constant  width  and  depth 
as  these  do  not  vary  in  any  systematic  way  along  the  channel,  and  includes 
a  region  of  tidal  flats  near  the  forcing  boundary.  The  area  of  flats 
decreases  linearly  to  zero  moving  downchannel.  Wachapreague  Inlet  (fig.  la) 
is  characterized  by  deep  (~5  m)  tidal  channels  which  terminate  in  broad 
shallow  bays.  The  bays  represent  extensive  tidal  flat  deposits  with  large 
changes  in  water  surface  area  over  a  tidal  cycle.  The  inlet  itself  is  rela¬ 
tively  stable  and  lacks  flood  tide  delta  deposits.  Wachapreague  is  modeled 
by  a  deep  channel  with  contracting  width  downchannel.  Channel  reaches  near 
the  interior  boundary  are  bordered  by  large  regions  of  tidal  flats.  Since 
these  experiments  do  not  represent  exact  matches  between  prototype  and  model 
geometry,  an  ensemble  of  runs,  varying  the  geometry  and  friction  parameters 
over  plausible  ranges  for  each  estuary,  was  considered  (table  8).  A  consis¬ 
tent  estuarine  response  over  a  reasonable  range  of  parameters  is  found  in 
each  case.  Two  particular  cases  (Nl,  W4)  were  selected  to  illustrate  the 
contrasting  features  of  estuarine  tidal  response  for  the  two  systems. 

The  models  representing  Wachapreague  are  characterized  by  a  longer  ris¬ 
ing  tide,  stronger  ebb  flows  and  relatively  weak  harmonic  growth  (fig.  29). 
This  result  holds  even  if  one  allows  for  a  trapezoidal  channel  geometry 
(table  8).  The  tidal  response  of  the  model  is  dominated  by  the  presence  of 
large  tidal  flat  area  and  a  relatively  small  a/h  ratio  ( ~0 .11).  Because 
of  these  factors,  even  strong  friction  (e.g.,  f  =  0.020)  is  unable  to  pro¬ 
duce  a  longer  falling  tide.  If  the  model  is  run  without  the  tidal  flats, 
harmonic  growth  is  weak  and  sea  surface  phase  produces  a  longer  falling  tide 
(table  8).  The  effect  of  the  tidal  flats,  as  indicated  by  the  previous 


TABLE  8 

NAUSET  MODEL  (L=7000  km) 
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Ratio  of  sea  surface  M4  and  M2  amplitudes  at  grid  point  U. 

Ratio  of  sea  surface  M4  and  M2  amplitudes  at  model's  interior  boundary. 
Relative  sea  surface  M2-M4  phase  at  model’s  interior  boundary. 

Relative  velocity  M2-M4  phase  at  model's  interior  boundary. 

Decline  in  tidal  range  at  model's  interior  boundary. 
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Figure  29.  Summary  of  the  tidal  response  of  Wachapreague  model  W4.  T 

panel  is  sea  surface  M2-M4  phase  adjustment.  Bottom  panel 
is  ratio  of  sea  surface  M4/M2  amplitudes. 
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section,  is  to  produce  both  greater  harmonic  growth  and  more  important,  a  M2 
M4  phase  in  sea  surface  (<j>c)  producing  a  longer  rise.  For  model  W4  rising 
tide  lasts  6  hours  25  minutes  and  falling  tide  6  hours  at  the  interior  bound 
ary.  Actual  field  measurements  at  Wachapreague  show  (M4/M2)<  equals  approxi 
mately  0.04  within  the  estuary  (Boon  and  Byrne,  1981).  Sea  surface  phase 
is  approximately  200-220°  producing  a  longer  rising  tide  and  stronger  ebb 
currents.  The  simple  model  described  in  this  section  reproduces  the 
general  features  of  this  tidal  response  by  focusing  on  the  "first  order" 
estuarine  characteristics  of  Wachapreague.  These  characteristics  include 
deep  tidal  channels  and  large  tidal  flats  at  the  distal  reaches  of  the 
channels.  Modeling  phases  and  amplitudes  exactly  would  require  (among 
other  factors)  a  representation  of  the  non-linear  hypsometric  curve  of  real 
estuaries  as  Boon  and  Byrne  (1981)  point  out.  The  model  presented  here 
employs  a  linear  hypsometric  curve. 

The  models  representing  Nauset  display  strong  harmonic  growth,  a  longer 
falling  tide  and  stronger  flood  flows  (fig.  30).  The  important  factors 
determining  its  tidal  response  are  the  large  a/h  ratio  (~0.30)  which  pro¬ 
duces  a  strongly  time-variable  channel  geometry  and  the  moderate  extent  of 
tidal  flats.  The  sea  surface  phase  relationship  (=64°)  indicates  that 
the  effects  of  channel  geometry  dominate  the  effects  of  tidal  flats.  A 
strong  time  asymmetry  exists  in  model  N1  with  falling  tide  lasting  7  hours 
15  minutes  and  rising  tide  5  hours  10  minutes  at  the  interior  boundary. 

Field  measurements  at  Nauset  (Chap.  3)  show  (M2/M4)/:  is  approximately  0.20 
and  <J>^,  63°.  Again,  a  simple  representation  of  estuarine  geometry  has 
reproduced  essential  features  of  Nauset 's  tidal  response.  The  important 
characteristics  in  this  example  are  shallow  tidal  channels  with  moderate 
tidal  flat  area  near  the  forcing  boundary.  The  tidal  distortion  in  sea 
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surface  within  the  channel  for  both  models  clearly  indicates  the  different 
time  asymmetries  developed  (fig.  31).  The  mean  velocity  signal  shows  the 
tendency  towards  ebb  dominant  flows  in  the  Wachapreague  case  and  flood 
dominant  in  the  Nauset  case  (fig.  32)  although  the  distortion  in  both 
instances  is  not.  exceedingly  large. 

5.  Discussion 

Although  the  modeling  is  a  simplified  representation  of  estuarine  tidal 
response,  the  results  can  be  combined  with  field  observations  of  asymmetry 
to  suggest  geometry  characteristics  producing  different  types  of  tidal  dis¬ 
tortion.  Harmonic  growth  in  both  models  and  natural  estuaries  is  a  strong 
function  of  the  non-linearity  parameter,  a/h.  Channel  geometry  effects 
become  important  only  when  this  parameter  is  relatively  large  ( >0 . 1 ) .  In 
this  case,  representing  channel  geometry  properly  (i.e.,  contracting  width 
with  depth)  has  an  important  effect  on  both  the  rate  of  harmonic  growth 
and  on  the  phasing  of  harmonics. 

Estuaries  characterized  by  relatively  large  values  of  th«  non-linearity 
parameter  (=0.3)  tend  to  have  longer  falling  tides  than  rising  tides  (flood 
dominant).  In  addition  to  Nauset,  Murrells  Inlet  which  has  approximately 
the  same  value  of  a/h,  also  has  a  longer  falling  than  rising  tide  (<t>:  =+90°). 
The  modeling  indicates  that  large  areas  of  tidal  flats  can  change  this 
behavior  to  a  situation  of  a  longer  rising  tide.  If  the  channel  geometry 
is  strongly  time-variable  (i.e.,  Ab/bo=0),  the  area  of  flats  must  be  more 
extensive  than  if  the  channel  were  rectangular.  It  is  difficult  to  put  a 
number  on  the  critical  area  of  tidal  flats  which  produces  a  longer  rising 
tide  because  the  modeling  is  not  an  exact  representation  of  estuarine 
physics  and  field  observations  are  sparse.  However,  most  of  the  very 
shallow  (a/h>~0.3)  estuaries  will  be  flood  dominant.  Deeper  estuaries 
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Figure  31a.  Distortion  of  surface  tide  at  grid  point  20  for  J  I  M  I 

model  runs  Nl  and  W4.  l\l  21  U  S  6  t  IT1  0  U  6  I  IMI 

- Wach.  model  W4 
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Periodogram  of  sea  surface  at  grid  point  20  for  model  runs 
N1  and  W4. 
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Figure  32.  Velocity  at  grid  point  6  for  model  runs  N1  and  W4.  -  N&US6t  fHO  (3  6 1  Nl 
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(a/h  ~0.1-0.2)  without  flats  will  develop  a  weak  tendency  towards  flood 
dominance.  The  presence  of  even  moderate  tidal  flat  area  can  reverse  this 
time  asymmetry.  Natural  estuaries  characterized  by  a  longer  rising  tide 
(e.g.,  North  Inlet,  SC;  Price  Inlet,  SC;  Wachapreague ,  V A)  have  relatively 
small  values  of  the  non-linearity  parameter  (~0.1)  and  extensive  regions 
of  tidal  flats.  Price  Inlet  and  Wachapreague  Inlet,  in  fact,  contain 
approximately  five  times  as  much  water  in  flats  as  in  channels  (Byrne  et 
al. ,  1975;  FitzGerald  and  Nummedal,  1983). 

Some  investigators  of  shallow  inlet/estuarine  systems  have  suggested 
that  the  variation  of  estuarine  surface  area  relative  to  channel  cross- 
sectional  area  controls  duration  asymmetries  (Nummedal  and  Humphries,  1978; 
FitzGerald  and  Nummedal,  1983).  This  argument  is  stated  utilizing  a  spati¬ 
ally  integrated  form  of  the  continuity  equation: 

dhb  =  Ac*U 
dt  Ab 

hb  =  sea  surface  level  in  estuary 
Ac  =  inlet  channel  cross-sectional  area 
Ab  =  area  of  estuary  at  hb 
U  =  mean  channel  velocity. 

The  spatial  integration  requires  estuarine  water  level  to  be  everywhere 
horizontal.  In  this  equation  the  rate  of  sea  level  change  is  proportional 
to  Ac/Ab.  The  ratio  is  a  maximum  at  low  tide  and  a  minimum  at  high  tide  in 
estuaries  with  large  changes  in  surface  area  (i.e.  extensive  tidal  flats). 
Consequently,  it  appears  that  the  water  surface  in  the  estuary  adjusts  more 
rapidly  to  the  ocean  tide  at  low  water  than  high  water.  The  channel  is 
said  to  be  more  "efficient"  in  its  exchange  of  water  at  low  than  high  tide. 
Essentially,  duration  asymmetries  are  explained  solely  by  continuity 


B.. 


-192- 


effects,  neglecting  momentum  considerations.  Ac/Ab  ratios  for  the  experi¬ 
ments  with  and  without  tidal  flats  exhibit  different  patterns:  channels 
without  flats  have  a  maximum  at  high  and  a  minimum  at  low  tide;  channels 
with  flats  show  the  opposite  trend  (table  9).  Nonetheless,  several  of  the 
latter  experiments  have  longer  falling  tides  than  rising  tides.  These 
experiments  are  characterized  by  strong  frictional  and  time-variable 
channel  geometry  effects  entering  both  momentum  and  continuity  equations. 

In  addition,  field  results  show  that  Nauset  Inlet  is  characterized  by  a 
longer  falling  tide  (Aubrey  and  Speer,  1983)  and  yet  Ac/Ab  is  a  minimum  at 
high  tide.  Clearly,  momentum  effects  cannot  be  ignored  in  all  cases. 
Shallow  frictionally  dominated  channels  develop  significant  pressure  gradi¬ 
ents  which  make  an  analysis  by  the  continuity  equation  incomplete.  This 
"continuity"  approach  is  best  suited  for  estuaries  which  consist  largely 
of  an  open  bay  that  oscillates  up  and  down  with  a  level  surface  in  response 
to  tidal  forcing. 

Although  the  modeling  focused  on  the  estuarine  response  to  M2  forcing, 
the  effects  of  harmonics  introduced  in  the  boundary  conditions  were 
studied.  The  question  is  how  important  is  estuarine  geometry  in 
determining  tidal  distortion  versus  harmonics  entering  the  estuary  through 
the  inlet.  Tidal  harmonics  can  be  generated  on  ebb-tide  deltas,  shallow 
platform  shoals  located  oceanward  of  the  inlet.  The  field  evidence  which 
exists  suggests  that  phase  relationships  and  amplitudes  of  harmonics  can 
be  strongly  modified  as  the  tide  propagates  into  the  estuary  (Boon  and 
Byrne,  1981;  Aubrey  and  Speer,  in  prep.).  At  Nauset,  pressure  sensor 
records  from  rhe  ebb-tide  delta  indicate  different  sea  surface  M2-M4  phase 
than  within  the  estuary  ($,-  ebb  delta  =30°;  tyr  estuary  =63°).  The  implica- 
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Exp  t.  13 

High  tide 
Mean  water 
Low  tide 


Expt.  11 

High  tide 
Mean  water 
Low  tide 


Expt.  17 

High  tide 
Mean  water 
Low  tide 


TABLE  9 
Ac/Ab  RATIOS 


level 


Ac/Ab  =  3.2  x  10"4 
Ac/Ab  =  2.5  x  10'4 
Ac/Ab  =  1.8  x  10~4 


level 


Ac/Ab  =  5.9  x  10‘5 
Ac/Ab  =  2.2  x  10~4 
Ac/Ab  =  2.7  x  10~4 


level 


Ac/Ab  =  1.3  x  10  4 
Ac/Ab  =  1.8  x  10"4 
Ac/Ab  =  1.8  x  10'4 
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tion  is  that  estuarine  geometry  is  a  strong  control  of  tidal  asymmetry,  and 
dissipation  within  the  inlet  proper  may  isolate  internal  harmonic  genera¬ 
tion  from  external  forcing. 

Two  different  values  of  4>c  (64°,  179°)  in  the  boundary  condition  and  the 
geometry  of  experiment  13  (table  2)  were  used  to  study  the  problem  in  a 
model  estuarine  channel.  One  of  the  phase  relationships  is  compatible 
with  the  channel  geometry  as  determined  by  the  runs  with  M2  forcing.  Intro¬ 
duction  of  a  small  amount  of  M4  (M4/M2  =  .033)  at  a  compatible  phase  in 
the  forcing  boundary  condition  eliminates  the  phase  adjustment  and  results 
in  slightly  greater  M4  growth  (fig.  33).  The  phase,  at  the  model's 
interior  boundary  is  close  to  that  of  experiment  13.  In  contrast,  the  same 
amount  of  M4  introduced  at  a  phase,  <j>,r,  of  179°  produces  a  different 
result.  The  M2-M4  phase  changes  dramatically,  approaching  the  value  of 
experiment  13  near  the  model's  interior  boundary,  and  the  amplitude  of  M4 
first  declines  before  growing  downchannel.  Thus  the  channel  geometry  in  the 
model  runs  exerts  an  important  control  on  tidal  asymmetry  regardless  of  the 
presence  of  harmonics  in  the  boundary  condition. 

The  implications  of  tidal  asymmetry  for  estuarine  sediment  transport 
are  clear.  Although  sediment  transport  represents  a  complex  flow/boundary 
interaction,  large-scale  patterns  of  estuarine  transport  are  controlled  by 
tidal  asymmetry.  Flood  dominant  estuaries  will  show  net  up-estuary  trans¬ 
port  if  the  asymmetry  is  sufficiently  pronounced.  If  sediment  supply  is 
adequate,  therefore,  these  systems  may  gradually  fill  over  time.  The  pro¬ 
cess  cf  shoaling  likely  will  accelerate  the  flood  dominant  tendencies  of  the 
of  the  estuary.  Eventually,  the  reduced  tidal  prism  from  infilling  can 
lead  to  inlet  shoaling  and  possibly  closure.  By  contrast,  deeper  systems 
with  tidal  flats  may  lepteseut  more  stable  estuarine  systems  with  ebb 
dominant  currents  able  to  flush  entering  sediments. 


W4  AMPLITUDE  (CM) 
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6.  Summary 

1) .  Numerical  integration  of  the  one-dimensional  equations  is  employed  to 

study  M2  tidal  propagation  in  long,  narrow  channels  typically  found  in 
shallow  estuaries  along  the  U.S.  east  coast.  Tidal  flats  are  fre¬ 
quently  found  along  the  estuarine  channels.  Numerical  results  indicate 
that  channels  without  tidal  flats  develop  a  time  asymmetry  character¬ 
ized  by  a  longer  falling  tide.  This  asymmetry  becomes  more  pronounced 
with  larger  friction  factors  and  increasingly  time-variable  channel 
cross-sectional  area.  Currents  associated  with  this  time  asymmetry 
tend  to  be  stronger  during  flood  than  ebb.  These  estuaries  can  be 
termed  flood  dominant. 

2) .  Addition  of  tidal  flats  to  these  channels  can  produce  the  case  of  a 

longer  rising  tide  and  stronger  ebb  currents  (ebb  dominant).  This 
time  asymmetry  develops  if  the  tidal  flat  area  is  large  enough  to  over¬ 
come  the  effects  of  time-variable  channel  geometry.  Channels  with 
relativey  small  values  of  a/h  and  rectangular  cross-sections  require 
smaller  tidal  flat  area  to  produce  an  ebb  dominant  asymmetry.  Weaker 
values  of  friction  with  tidal  flats  can  also  drive  the  time  asymmetry 
towards  a  longer  rising  tide. 

3) .  The  general  features  of  tidal  channels  in  shallow  estuaries  can  be 

represented  with  simple  geometries.  Despite  the  physical  complexity 
of  these  systems,  essential  features  of  the  estuarine  tidal  response 
in  sea  surface  can  be  recovered  from  simple  models.  This  is  illus¬ 
trated  by  modeling  idealized  representations  of  tidal  channels  at 
Nauset  Inlet,  MA  and  Wachapreague  Inlet,  VA.  Over  a  range  of  geometry 
parameters  and  friction  factors  which  could  characterize  these  estuar¬ 
ies,  reasonable  agreement  in  M4/M2  amplitude  ratios  and  M2-M4  phas¬ 
ing  in  sea  surface  is  demonstrated. 
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h) .  Field  observations  and  model  results  suggest  that  the  shallowest 

estuaries  (a/h>0.3)  generally  will  be  characterized  by  longer  falling 
tides  and  stronger  flood  currents  unless  extensive  tidal  flats  are 
present.  Deeper  systems  (a/h~0 . 1- 0 . 2 )  with  tidal  flats  will  likely 
have  longer  rising  tides  and  stronger  ebb  currents.  Sediment 
transport  patterns  in  flood  dominant  estuaries  will  be  directed 
up-estuary  with  consequent  channel  shoaling  if  sediment  supply  is 
adequate.  Ebb  dominant  estuaries  may  represent  more  stable 
conf igurations  with  stronger  ebb  currents  able  to  flush  sediments 
entering  the  estuary  through  the  inlet. 

5).  The  models  considered  here  provide  only  qualitative  information  on  the 
effects  of  tidal  flats.  To  include  these  effects  more  accurately,  the 
modeling  must  be  extended  to  two  dimensions  (depth-integrated).  Two- 
dimensional  models  solve  the  full  momentum  and  continuity  equations 
over  the  flats  as  well  as  in  the  channel.  This  enables  the  impact  of 
flats  in  momentum  to  be  more  carefully  evaluated.  Allowing  greater 
cross-sectional  area  variability  along  channel  and  the  possibility  of 
channels  in  the  tidal  flats  also  requires  a  two-dimensional  model. 
Given  sufficient  spatial  resolution  of  estuarine  geometry,  two-dimen¬ 
sional  models  can  also  more  accurately  include  regions  where  the  advec 
tive  term  is  of  importance  (e.g.,  channel  bends,  channel  constric¬ 
tions).  Future  efforts  wilL  be  directed  towards  this  type  of  modeling 
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Ar PEND IX  I 

EQUATIONS  OF  MOTION  FOR  ONE-DIMENSIONAL  FLOW 
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1.1.  Introduction 

The  channel  geometry  considered  in  this  derivation  is  a  narrow  and 
shallow  channel  with  a  well-mixed  water  column  (fig.  1.1).  The  axes  are 
as  shown  in  the  figure  and  the  geometry  suggests  that  the  longitudinal 
velocity  component,  u,  will  v -  much  greater  than  the  lateral,  v,  and 
vertical,  w.  The  derivation  assumes  that  changes  in  width  and  cross- 
secLional  area  are  gradual  in  the  longitudinal  direction.  The  dynamics 
will  be  expressed  in  terms  of  the  volume  flux,  U(x,t),  cross-sectionally 
averaged  velocity,  u(x,t),  and  surface  elevation,  q(x,t).  The  derivation 
presented  below  follows  Dronkers  (1964),  Madsen  (unpub.  MIT  class  notes) 
and  Brown  and  Trask  (1980). 

1.2.  Continuity  Equation 

The  continuity  equation  for  an  incompressible  fluid  is: 

3u  +  3 v  +  3w  =  0  (1) 

3x  3y  3z 

This  equation  is  integrated  over  the  cross-section  shown  in  Figure  1.1. 


b2  C 

s  x 

-bi-h 


3u  +  3w 
3x  3y 


dzdy  =  0 


(2) 


b2C  b2  C 

In  the  following,  the  integrals  X  X  dzdy  ,  X  dy  and  X  are  denoted  <...>, 

-bi-h  -bi  -h 

{...}  and  ...  respectively.  Evaluation  of  the  integrals  in  (2)  requires 


repeated  use  of  Liebnitz'  rule: 


<3u>  =3  <u>  -  3b2  u2  -3bi  ui  +  {u(-h)3(-h)  -u(£)9£}  (3a) 


3x 

3x 

3x  3x 

3x 

3x 

<3v> 

3y 

=  V2  -  Vi 

+  {v(-h)3(:ih) 
3y 

-  v(U9i) 

3y 

(3b) 

<  3w  > 
3z 

=  {w(U  - 

w(-h ) } 

(3c) 

PLAN  VIEW 

Figure  1.1.  Channel  geometry  and  coordinate  system. 
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The  subscripts  1,2  refer  to  conditions  at  y=-bi,b2  respectively.  To 
simplify  these  integrals,  kinematic  surface  and  bottom  boundary  conditions 
are  applied: 


D£  =  w(C) 

(4a) 

Dt 

D(-h)  =  w(-h) 

(4b) 

Dt 

and  the  integrated  continuity  equation  becomes: 

3  <u>  +  {3£  -  _3(-h)}  +  v 2  -3b2v2  -  vi  -  3b!V2  =  0  (5) 

3x  3t  3t  3x  3x 

The  integral  {3£  -  3(-h)}  can  be  rewritten  using  the  Liebnitz  rule  as, 
3t  3t 

_3{C  -  (-h)}  -  3b2(S  -  (-h))2  -  3bi (C  -  (~h) ) i  (6) 

3t  3t  3t 

and  utilizing  kinematic  side  boundary  conditions, 

Dbz  *  v(b2)  (7a) 

Dt 

Dbi  =  v(bi )  (7b) 

Dt 

the  integrated  continuity  equation  can  be  written: 

a  <u>  +  a  {c  -  (-h)}  =  0  (8) 

3x  at 

Defining, 

<u>  ~  U  =  u*A  A  =  cross-sectional  area 

{£  -  (-h)}  =  A 
continuity  becomes 

3  U  +  3A  =  0  (9) 

3x  at 

For  the  trapezoidal  cross-section  shown  in  figure  1.2, 

3A  =  b3£  b  =  channel  width  (10) 

3t  3t 
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The  forms  of  continuity  considered  in  the  numerical  model  result  from 
substituting  (10)  into  (9): 


9U  + 

b9£  =  0 

(11a) 

9x 

9t 

9A*u 

o 

ii 

SI 

+ 

(lib) 

9x 

9t 

Note  that  the  derivation  does  not  account  for  lateral  inflow  from  side 
channels . 

I. 3  Momentum  Equation 

The  geometry  of  the  problem,  asumptions  of  gradual  change  in  cross- 
section  in  the  longitudinal  direction  and  hydrostatic  pressure  distribution 
allow  us  to  state  the  y-  and  z-momentum  equations  as: 

C  =  C(x,t)  (12) 

p  =  pg(C-z)  (13) 

The  x-momentum  equation  is 

8u  +  u9u  +•  w9u  =  -  1  9p  +  1  Qxx  x  +  9xvx  +  9t vz  )  (14) 

9t  9x  9y  p  9x  p  9x  9y  8z 

Adding  to  (14),  u(V  •  v),  where  v(u,v,w)  is  the  velocity  vector,  and  substituting 
(13)  in  (14),  x-momentum  becomes 

9u  +  9u2  +  9uw  =  -g9£  +  l  (9x x x  +  9xvx  +  9xvz  11 .  (15) 

9t  9x  9y  9x  p  9x  8y  9z 

The  velocities  in  (15)  are  time  averaged  over  an  interval  long  compared  to 
turbulent  time  scales  but  short  compared  to  tidal  time  scales.  The  shear 
stresses,  x,  represent  both  viscous  forces  and  stresses  arising  from  advec- 
tion  of  momentum  by  turbulent  velocity  fluctuations  (Reynold's  stresses). 

It  has  been  assumed  in  deriving  (15)  that  the  coriolis  term  is  not 
important  in  a  narrow  tidal  channel.  This  assumption  is  valid  if  the  Kelvin 
number,  Ke  =  fb  (gh)~1/2,  is  much  smaller  than  unity.  For  the  estuarine 
system  of  interest,  b~200  m,  h~4  m,  f-lO’^s,  and  KE~10"3. 
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The  momentum  equation  (15)  is  integrated  over  the  channel  cross-section. 

This  integration  yields  for  the  total  derivative  irr, 

<Du>  =  _8  <u>  +  _8  <u2>  -  {u(C)(3£  +  u(C)„8£  +  v(£)8£  ~  w(C)} 

Dt  3t  3x  8t  3x  8y 

+  {u(-h)(8(:^h)  +  u(-h)8(=h)  +  v(-h)8(::h)  -  w(-h)} 

8t  8x  8y 

c  c 

-J"  U2(8b2  +  U28t^2  -  v2)dz  -  f  ui(8b_^  +  ui8b_^  -  V[)dz.  (16) 

-h  8t  8x  -h  8t  8x 

The  third  and  fourth  integrals  vanish  by  applying  kinematic  boundary  condi¬ 
tions  on  the  free  surface  and  on  the  channel  bottom.  If,  in  addition,  no 
lateral  inflow  takes  place,  the  last  two  integrals  vanish  due  to  kinematic 
conditions  on  the  channel  side.  Therefore, 

<Du>  =  _8  <u>  +  _8  <uz>  (17) 

Dt  8t  8x 

The  integral  _8  <uz>  can  be  written  as: 

8x 

_8  <uz>  =  _8  Km  U2  (18) 

8x  8x  A 

Km  =  momentum  coefficient  =  <uz>  ~  0(1) 

U*u 

The  momentum  coefficient  accounts  for  the  deviation  from  constant  velocity 
over  the  cross-section  and  allows  the  advective  term  to  be  written  directly 
in  terms  of  the  volume  flux,  U.  For  gentle  variation  in  cross-sectional 
area.  Km  is  taken  equal  to  unity  and. 


<  u 2  >  =  u2 
A 

The  acceleration  term  becomes,  therefore, 

<Du>  =  8U  +  _8  U2 
Dt  8t  8x  A 


(19) 
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This  derivative  can  be  expressed  in  terms  of  the  mean  velocity,  u,  by  manipul¬ 
ating  the  definitions  of  U  and  A: 

9  <u>  =  A  •  9u  +  A  •  u9u  -  9_(1<u2>)  (20) 

9t  3t  9x  3x  A 

Dividing  through  by  A  and  substituting  in  (17)  yields: 

<Du>  =  9u  +  u9u  +  1  9_(<u2>  -  <u>2)  (21) 

Dt  3t  3x  A  3x  A 

The  third  term  in  (21)  accounts  for  the  deviation  from  constant  velocity 
over  the  cross-section  and  will  be  small  under  the  same  conditions  that  Km 
approximates  unity. 

Assuming  a  gentle  lateral  and  depth-wise  variation  over  the  cross- 

section,  the  pressure  term  in  (15)  can  be  integrated: 

-<g9£>  =  -gA9£  (22) 

9x  3x 

The  shear  stress  terms  are  integrated  after  neglecting  the  term  9txx . 

3x 

This  requires  the  assumption  of  slow  longitudinal  variations  in  channel 
geometry.  Integration  of  the  remaining  terms  yields 

(1( 9t x v  +  9t x z ) ^  =  Xxy»  2  —  Txy,  i  +  {tXz (( )}  — 

p  3y  3z 

{tXI(-h)  -  9_(-h)xxy(-h)}  (23) 

3y 

The  first  two  integrals  represent  sidewall  friction  if  channel  sidewalls 
are  vertical.  The  channel  geometry  under  consideration  (fig.  1.1)  does 
not  have  steep  lateral  boundaries  and  the  sidewalls  will  be  considered 
part  of  the  channel  bottom  boundary.  Thus  the  first  two  integrals  are 
neglected.  For  the  case  of  a  rectangular  channel,  the  lateral  stresses 
which  scale  as  h  are  unimportant  compared  to  the  bottom  stresses  which 


scale  as  b. 
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The  third  integral  represents  a  surface  wind  stress  and  will  be  taken 

as  zero  for  the  tidal  problems  under  consideration.  The  final  integral 

can  be  shown  to  be  equivalent  to  an  average  shear  stress  on  the  solid 

boundary  times  the  wetted  perimeter: 

{t:« (-h)  -  3_(-h)Txy(-h)}  =  -  Tb*P  (24) 

3y 

Tb  =  average  shear  stress 
P  =  wetted  perimeter. 

The  momentum  equation  is  closed  mathematically  by  writing  shear  stress  in 
terms  of  mean  velocity,  u,  or  volume  flux,  U: 


Tb  =  pf 

|u|u 

(25a) 

Tb  =  pf. 

|u|u 

(25b) 

A«A 

Combining  (25)  and  (22)  with  (21)  or  (19)  yields  the  final  form  of  the 
momentum  equation: 


3U 

+  3_  U2 

=  -gA3£  -  f | U | U  •  P 

(26a) 

3t 

3x  A 

3x 

A*A 

3u 

+  u3u  = 

-g9£  -  f lu|u  •  P 

(26b) 

3t 

3x 

3x 

A 

1.4  Tidal  Flats 

The  effect  of  lateral  inflow  from  tidal  flats  on  the  continuity  and 
momentum  equations  can  be  included  in  a  more  systematic  manner  than  shown 
in  Chapter  4.  If  lateral  inflow  into  the  main  channel  occurs,  the  kine¬ 
matic  side  boundary  conditions  (7a, b)  do  not  vanish  and  the  continuity 
equation  (11a)  is  rewritten: 

3U  +  b8£  +  qE  =  0  (27) 

9x  3t 

qe  =  lateral  inflow  per  unit  length  of  channel 
b  =  width  of  momentum  transporting  channel 
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The  amount  of  inflow  (or  outflow)  caused  by  tidal  changes  in  sea  surface  is 
given  by: 

qE  =  (bT  -  b)3£  (28) 

3t 

bT  =  total  width  of  channel  including  flats 
substituting  (28)  into  (27),  the  equation  of  continuity  becomes 

au  +  bTay  (29) 

3x  at 

which  is  the  same  as  equation  (6)  in  Chapter  4. 

Similarly  for  momentum,  the  retention  of  lateral  inflow  terms  in  (26b) 
allows  the  equation  to  be  written 

3u  +  u3u  +(ue  -  u)qE  =  -g3£  -  f |u|u  •  P  (30) 

3t  3x  3x  A 

ue  =  longitudinal  velocity  component  of  lateral  inflow 

If  the  lateral  inflow  from  tidal  flats  has  approximately  zero  velocity, 

ue  -  0,  and  substituting  (28)  and  (29)  into  (30)  yields: 

3u  +  u3u  +  u(bt  -  b)3U  =  -g3£  -  f  juju  •  P  (31) 

3t  3x  bT  3x  0x  A 

In  terms  of  discharge,  U,  momentum  is  written 

au  +  au2  +  (br  -  b0 )U  au  =  -gA3£  -  f lulu  •  P  (32) 

3t  3x  A  bT  A  3x  3x  A«A 

where  the  lateral  inflow  term  can  be  compared  to  equation  (10)  in  Chapter  4. 
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